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1 6 Abstract 


The literature has shown that crack propagation in cracked metal sheets can 
be significantly reduced by bonding an uncracked reinforcement to the metal 
sheet. However, cyclic debonding typically occurs over a localized area near the 
crack. Herein, an analysis was developed to predict both the crack growth and 
debond growth in a reinforced system. The analysis was based on the use of com- 
plex variable Green's functions for cracked, isotropic sheets and uncracked, 
orthotropic sheets to calculate Inplane and interlaminar stresses, stress inten- 
sities, and strain-energy-release rates. An iterative solution was developed 
that used the stress intensities and strain-energy-release rates to predict crack 
and debond growths, respectively, on a cycle-by-cycle basis. The analysis was 
verified with experiments. 

The analysis was used in a parametric study of the effects of boron-epoxy 
composite reinforcement on crack propagation in aluminum sheets. The study 
showed that the size of the debond area has a significant effect on the crack 
propagation in the aluminum. For small debond areas the crack propagation rate 
is reduced significantly, but these small debonds have a strong tendency to en- 
large. Debond growth is most likely to occur in reinforced systems that have a 
cracked metal sheet reinforced with a relatively thin composite sheet. 


17 


The analysis predicts crack growth in reinforced systems. Hence, the analy- 


sis can be applied in developing methods 
Increase the lives and payloads of metal 

to repair damaged metal structures and to 
structures by selective reinforcement. 
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CRACK PROPAGATION IN ALUMINUM SHEETS 


REINFORCED WITH BORON-EPOXY 
by 

G. L. Roderick 
ABSTRACT 

The literature has shown that crack propagation in cracked 
metal sheets can be significantly reduced by bonding an uncracked 
reinforcement to the metal sheet. However, cyclic debonding typically 
occurs over a localized area near the crack. Herein, an analysis was 
developed to predict both the crack growth and debond growth in a 
reinforced system. The analysis was based on the use of complex 
variable Green's functions for cracked, isotropic sheets and uncracked, 
orthotropic sheets to calculate inplane and interlaminar stresses, 
stress intensities and strain-energy-release rates. An iterative 
solution was developed that used the stress intensities and strain- 
energy- release rates to predict crack and debond growths, respectively, 
on a cycle-by-cycle basis. The analysis was verified with experiments. 

The analysis was used in a parametric study of the effects of 
boron-epoxy composite reinforcement on crack propagation in aluminum 
sheets. The study showed that the size of the debond area has a 
significant effect on the crack propagation in the aluminum. For 
small debond areas the crack propagation rate is reduced significantly, 
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but these small debonds have a strong tendency to enlarge. Debond 
growth is most likely to occur in reinforced systems that have a 
cracked metal sheet reinforced with a relatively thin composite sheet. 

The analysis predicts crack growth in reinforced systems. 

Hence, the analysis can be applied in developing methods to repair 
damaged metal structures and to increase the lives and payloads of metal 
structures by selective reinforcement. 



TABLE OF CONTENTS 


ABSTRACT i 

TABLE OF CONTENTS iii 

LIST OF FIGURES vi 

LIST OF TABLES ix 

LIST OF SYMBOLS x 

CONVERSIONS FROM CUSTOMARY TO SI UNITS xvii 

Chapter 

I. THE USE OF COMPOSITE REINFORCEMENT TO PREVENT FATIGUE 

FAILURE IN AIRCRAFT 1 

II. FATIGUE BEHAVIOR OF THE REINFORCED SYSTEM CONSTITUENTS . . 5 

Fatigue of Metals 5 

Fatigue of Bonded Systems 10 

Fatigue of Composite Materials 12 

III. FATIGUE TESTS OF THE REINFORCED SYSTEM 16 

IV. STRESS ANALYSIS 24 

Formulation of Linear Elastic Solution 24 

Numerical Solution 37 

Nonlinear Solution 41 

V. FATIGUE ANALYSIS 49 

Crack Growth Rate 49 

Debond Growth Rate 52 

Prediction of Crack and Debond Growth 60 

iii 



VI. ASSESSMENT OF ACCURACY 62 

Numerical Integration 62 

Accuracy of the Analysis 71 

VII. PARAMETRIC STUDIES ON CRACK AND DEBOND GROWTH 79 

Stress Intensity Factor, Crack Growth Rate 84 

Strain Energy Release Rate, Debond Propagation 84 

Nonlinear Effects 91 

Prediction of Crack and Debond Growth 94 


CONCLUSIONS 95 

APPENDICES 

A. DETERMINATION OF DEBOND CONSTANTS 97 

Specimen Fabrication 97 

Fatigue Tests 100 

Stress Analysis 102 

Calculation of Strain Energy Release Rates 108 

Curve Fit for Empirical Constants 112 

B. ADHESIVE SHEAR DEFORMATION ASSUMPTION 116 

Adhesive Bulk Properties 116 

One-Dimensional Solution 121 

Finite Element Solution 126 

One-Dimensional Versus Finite Element Solution 128 

C. REMOTE STRESSES IN THE ADHERENDS 134 

D. GREEN'S FUNCTIONS FOR THE CRACKED SHEET 137 

Green's Functions for Stress 161 

Green's Functions for Displacements 166 


IV 



E. GREEN'S FUNCTIONS FOR AN UNCRACKEO ORTHOTROPIC SHEET. . . 173 

Green's Functionsfor Stresses 179 

Green's Functionsfor Displacements 181 

F. COMPUTER PROGRAM 185 

Sample Input 189 

Sample Output 189 

Program Listing 193 

LIST OF REFERENCES 224 


V 



LIST OF FIGURES 


Figure Page 

1. Stress Distributions near Crack Tip 7 

2. Panel Configuration 17 

3. Crack Growth Versus Applied Load Cycles 20 

4. Experimental Crack Propagation Rates 21 

5. Ultrasonic C-scans of Fatigued Specimens 22 

6. System to be Analyzed 25 

7. Inplane and Interlaminar Stress in the Reinforced 

System 27 

8. Domain of Integration 38 

9. Freebody for Determination of G along Debond Front . 53 

10. Contour for Determination of G for Strip 55 

11. Flow Chart for Fatigue Analysis 61 

12. Domain B for Convergence Analysis 65 

13. Convergence of Interlaminar Shear Stresses 67 

14. The Effect of Mesh Size on Interlaminar Stresses ... 70 

15. Mesh for Analysis of Panels A and B 73 

16. Experimental and Calculated Crack Lengths as a Function 

of Applied Load Cycles 74 

17. Experimental and Calculated Crack Propagation rates. . 76 

18. Debond Aspect Ratio Versus Crack Length 77 

19. Stress Intensities for 0.05 in Thick Metal Adherend. . 85 

vi 



Figure Page 

20. Stress Intensities for 0.10 in Thick Metal Adherend. . 86 

21. Stress Intensities for 0.15 in Thick Metal Adherend. . 87 

22. Strain Energy Release Rate for 0.05 in Thick Metal 

Adherend 88 

23. Strain Energy Release Rate for 0.10 in Thick Metal 

Adherend 89 

24. Strain Energy Release Rate for 0.15 in Thick Metal 

Adherend 90 

A 1. Debond Specimen Configuration 98 

A. 2. Debonding Behavior 101 

A 3. Beam Model for Stress Analysis 104 

A. 4. Comparison of Finite Difference Solution with 

Experimental Data 107 

A. 5. Strain Energy Release Rate Versus Debond Front 

Location Ill 

A. 6. Correlation of Debonding Rates with Strain Energy 

Release Rates 114 

B. l. Bulk Adhesive Specimen Configuration and 

Instrumentation 117 

3.2. Adhesive Stress Strain Curve 119 

B.3. Specimen Configuration and Freebody For One- 

Dimensional Solution 122 

B.4. Finite Element Mesh 127 

B.5. Finite Element Versus One-Dimensional Solution .... 129 

D.l. Superposition Method used to Formulate Green's 

Functions for a Cracked Sheet 139 

D.2. Path for Contour Integration 146 

D. 3. Finite Element Mesh used to Check Green's Functions . 165 


VI 1 



Figure Page 

D.4. Verification of Green's Functions for Stresses in a 

Cracked Isotropic Sheet 167 

D. 5. Verification of Green's Functions for Displacements 

in a Cracked Isotropic Sheet 172 

E. l. Locations of Point Loads on an Orthotropic Sheet . . . 178 

E.2. Verification of Green's Functions for Stresses in an 

Uncracked Orthotropic Sheet 182 

E. 3. Verification of Green's Functions for Displacements 

in an Uncracked Orthotropic Sheet 184 

F. l. Model for Sample Run 188 


viii 



LIST OF TABLES 


Table Page 

1. CRACK LENGTHS AND PROPAGATION RATES 19 

2. CALCULATED VALUES FOR PARAMETRIC STUDIES 80 

A. l. STRAIN ENERGY RELEASE RATES AND DEBOND PROPAGATION RATES . 113 

B. l. EFFECTIVE SHEAR MODULI FOR REINFORCED SYSTEM 133 


IX 



LIST OF SYMBOLS 


a 


da/dN 

A(C) 





1 1 


b 

db/dN 

bb 

BB 


c 


I 


Half-length of crack, in 

Crack propagation rate in metal, in/cycle 

Kernel for contour integration 

Constant coefficients in stress functions for 
orthotropic sheet 

Expressions used for determination of yielding of the 
adhesive layer 

Constants for analysis of unsymmetrical beams under 
axial loads 

Integrals of displacement Green's functions for cracked 
metal sheet, in 

Debond length along the y-axis, in 

Debond propagation rate along the y-axis, in/cycle 

Integrals used to find displacements from interlaminar 
shear stresses 

Integrals of displacement Green's functions for 
orthotropic sheet, in 

Complex function used in displacement Green's functions 
for cracked metal sheet, in 

Constant coefficient in crack propagation equation for 
cracked metal sheet, in 


X 



c 

2 

[C] 

cn, C12 

C21, C22 
d 

[D] 


e 

E 

f 

F 


F 


P 


F , 
1 


F 

2 


F , 

3 


F 

4 


F . F 

S 6 


Constant coefficient in debond propagation equation 
Stiffness matrix for an isotropic sheet in plane stress, 
psi 

Complex constants used in orthotropic stress 
analysis 

Distance from test machine grips to change in cross 
section of debond specimen 
Stiffness matrix for an orthotropic sheet in plane 
stress, psi 

Numerical constant 2.718281824**** 

Eccentricity of axial load in beam column analysis, in 
Young's modulus of elasticity, psi 
Interlaminar shearing stresses, psi 
Aspect ratio of elliptical debond b/a 
Complex function that relates pressure on a crack to 
displacement field, psi 
Complex function that relates point loads to 

displacements in an uncracked isotropic sheet, in 
Complex functions that relate remote stresses in a 
cracked metal sheet to displacements, in 
Complex functions that relate remote stresses to 
displacements in an uncracked orthotropic sheet, in 
Complex functions that relate interlaminar shear 
stresses to displacements in a cracked isotropic sheet, 
in 


xi 



g 


9o 


G 

G 

^eff 

Gl, G1A, G2, 


G3, G5, G6, G7 
G8, G9 


GD 


GS 


HD 

HS 

Kz) 

II(z) 


Complex functions that relate interlaminar sheet 
stresses to displacements in an uncracked orthotropic 
sheet, in 

Increment of interlaminar stresses between loads that 
cause successive yielding of elements, psi 

Complex function used in developing the displacement 
Green's functions for a cracked isotropic sheet 

Strain energy release rate, in-lb/in 

Shear modulus, psi 

Effective shear modulus of the adhesive, psi 

Complex functions used to express stress Green's 
functions for cracked isotropic sheet, psi 

Complex stress functions used to express stress Green's 
functions for an uncracked orthotropic sheet, psi 

Displacement Green's functions for a cracked isotropic 
sheet, in 

Stress Green's functions for a cracked isotropic 
sheet, psi 

Displacement Green's functions for an uncracked 
orthotropic sheet, in 

Stress Green's function for an uncracked orthotropic 
sheet, psi 

Complex root of /z^ - eP” defined by ,, 

1 im I (z) 

z -> ® ~z ^ ’ 

Complex integral used in developing displacement Green's 


functions for a cracked isotropic sheet 



Path independent integral used to determine the 
strain energy release rate 
Constant used in Von Mises yield criterion 
Constant for one-dimensional analysis 
Stress intensity for mode I, psi-in 
Stress intensity for mode II, psi-in^ 

u 

Critical mode I stress intensity at failure, psi-in 
Stress intensity factor 

Ratio of the shear modulus and Young's modulus for the 

orthotropic sheet G /E 

xy y 

Moment in beam column analysis, in- lb 
Exponent for crack propagation equation for metal sheet 
Exponent for debond propagation equation 
Ratio of the Young's moduli of the orthotropic sheet, 
Ex/Ey 

Complex constants for orthotropic stress analysis 
Normal stress acting on crack surface, psi 
Axial load used in beam column analysis, lbs 
Complex constants for orthotropic analysis 
Shear stress acting on crack surface, psi 
Polar coordinates, in, deg 
Small radius for contour integration, in 
Ratio of minimum to maximum load in a load cycle 
Large radius for contour integration, in 
Remote stress applied to the reinforced system in the 
y-direction, psi 



S , s 
1 2 


S 


t 


u, V, w 


U 

V 


w , w 
1 2 


w 


W 


e 


X, y, z 


X 


0 ’ 


^0 


X, Y 


XA, XB, 
XC, XI 
XK 


z 

Z , z 
1 2 

a 

6 

y 


Complex roots from the governing equation for 
orthotropic sheets in plane stress or strain 
Complex representation of point load acting in complex 
plane per sheet thickness Ibs/in 
Thickness of adherends or adhesive, in 
Displacements in cartesian coordinates in the x, y, and 
z-directions respectively, in 
Internal strain energy, in-lbs 
Shear in beam, lbs 

Complex variables for orthotropic analysis 
External work, in- lbs 
Strain energy density 
Cartesian coordinates 

Cartesian coordinates of point of load application 
in the complex plane, in 
Unit load components in the x and y-directions 
respectively, lbs 

Complex functions used in the displacement Green's 
functions for a cracked isotropic sheet, in 
Complex function used as Green's function for stress 
intensities, psi-in^ 

Complex variable for analysis of isotropic sheet 
Complex variables for analysis of orthotropic sheet 
Constant for shear lag analysis, in 
Constant for shear lag analysis 
Shear strain in the adhesive 


XIV 



6 

A 

e 

C 

00 

V 

IT 

a 

T 

4 > 

$ , $ 

1 2 

U3 

Q. 

n 

Subscripts 

ad 

c 

m 

xy 


Axial deformation, in 

Operator to denote small change in quantities 
Strain 

Complex variable for contour integration 
Rotation in beam 
Poisson's ratio 
Constant 3-1415926535 
Stress, psi 
Stress in metal, psi 
Stress in composite, psi 
Shear stress in the adhesive, psi 
Stress function for isotropic sheet 
Derivative of stress functions for isotropic sheet, (j' 
Stress functions for orthotropic sheet 
Stress function for isotropic sheet 
Derivative of stress function for isotropic sheet, ij;' 
Stress function for isotropic sheet 
Derivative of stress function for isotropic sheet, w' 
Constant for plane stress or strain used in complex 
variable formulation 

Adhesive 

Composite adherend 
Metal adherend 

Quantity in the x-direction due to a load in the 


y-direction 


XV 



XX Quantity in the x-direction due to a load in the x- 

direction 

yx Quantity in the y-direction due to a load in the 

x-di recti on 

yy Quantity in the y-direction due to a load in the 

y-direction 


Superscripts 

Conjugate value, if z = x + iy then z = x - iy 
and if F(z) = A + iB then F(z) = A - iB and 
F(z) = F(I)'. 


XVI 



CONVERSIONS FROM CUSTOMARY TO SI UNITS 


FORCE 

Pound force (1b) = 4.448 Newtons (N) 

PRESSURE 

Stress (psi) = 6894 Pascals (Pa) 

DISTANCE 

Inches (in.) = 0.0254 meters (m) 

ENERGY 

Foot-pound force (ft-lb) = 1.355 Joules (J) 
TEMPERATURE 

Fahrenheit (°F) = 5/9 (tp - 32) Celsius (°C) 



CHAPTER I 


THE USE OF COMPOSITE REINFORCEMENT 
TO PREVENT FATIGUE FAILURE IN AIRCRAFT 

A potential cause of aircraft crashes is fatigue failure. As 
shown by Hardrath (1971), most types of civil aircraft have experienced 
some form of fatigue problem. In addition Lowndes and Miller (1969) 
indicate that fatigue failures have frequently occurred in military 
aircraft. In some cases the fatigue failures led to loss of lives and 
the aircraft. In efforts to eliminate such failures, both government 
research laboratories and aircraft manufacturers have studied the 
fatigue failure process in depth. These studies showed that the rate 
at which the fatigue damage develops in metals is a function of the 
stress level in the structure and occurs in three stages: crack 

initiation, stable crack propagation, and unstable crack propagation 
(catastrophic failure). Although aircraft structures can be designed 
to have stresses low enough to prevent fatigue failures, the weight 
penalty would be enormous and would make the aircraft uneconomical 
to operate. Hence, a trade-off exists between low stresses and low 
weight, and weight efficient structures will almost always have 
stresses high enough to support fatigue damage accumulation. 

Fatigue cracks initiate at local stress concentrations in the 
structure The local stress concentrations may be caused by poor 
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fatigue design, by manufacturing defects in the material, or by damage 
caused by the flight environment. Although methods can be employed to 
reduce the occurrence of fatigue crack initiation, the development of 
such cracks seems almost inevitable. 

Once a crack initiates it grows at a stable rate until it reaches 
some predictable, critical length after which catastrophic failure follows. 
Fortunately, in aluminum aircraft structures the critical crack length 
is large and the crack is easy to detect long before it reaches a 
critical length. Consequently, fatigue cracks can be tolerated in 
an aircraft structure as long as the structure is inspected periodically 
to locate cracks before they become critical. Of course, once the 
crack is detected it must be repaired before it becomes critical. The 
repair can be made by either replacing the component or by repairing 
it in situ. Because replacing a component may involve high cost and 
keep the aircraft out of service for extended periods of time, repair- 
ing the component in situ is frequently very desirable. 

Basically, a fatigue crack can be repaired by reducing the stress 
state in the vicinity of the crack tip. One method of reducing the 
stress state is to reinforce the crack with unidirectional composite 
(fibers are perpendicular to plane of crack). The composite reinforce- 
ment reduces the stress state near the crack tip by two mechanisms. 

First, adding the composite reinforcement lowers the overall stress in 
the cracked metal by increasing the cross-sectional area and by 
providing an alternate, stiffer load path by virtue of the high 
modulus of the composite. This reduction in stress can be easily 
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calculated by simple strength of materials theory and, hence, is 
easy to investigate. The second mechanism comes from the development 
of stresses between the metal and composite adherends. Several papers 
(Kula et al 1973, Ellis 1976, Johnston and Stratton 1975, Ratwani 1977) 
have shown that these interlaminar stresses have a profound effect on 
crack propagation. These interlaminar stresses reduce the stresses 
near the crack tip and, consequently, retard its growth. As will be 
shown later, the investigation of these interlaminar stresses requires 
a much more extensive analysis than that provided by strength of 
materials theory. 

Consequently, a need exists for the development of a realistic 
fatigue analysis that incorporates the effects of the interlaminar 
stresses. Accordingly, the objective of this dissertation was to 
develop such an analysis and use it to study the effects of composite 
reinforcement on the fatigue life of cracked metallic structure. 

To meet this objective the following approach was taken. First, 
in Chapter II the fatigue behavior of the constituents of the reinforced 
system was characterized. Next in Chapter III the fatigue behavior of 
the reinforced system was studied experimentally. Then, in Chapter IV 
with the use of the results of Chapters II and III and complex variable 
theory, a static analysis was developed that related applied loads, 
adherend thicknesses, debond size, and crack length to crack propagation 
rates. Next, in Chapter V the analysis was further developed to predict 
both debond and crack growth as a function of applied load cycles. The 
accuracy of the analysis was investigated in Chapter VI. Finally, in 
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Chapter VII the analysis was used to parametrically study crack 
propagation in reinforced systems. 

Throughout the development of the analysis many items required 
detailed analytical or experimental investigation. These investigations 
were developed and discussed in several appendices. 



CHAPTER II 


FATIGUE BEHAVIOR OF THE REINFORCED SYSTEM CONSTITUENTS 

The reinforced system that will be considered herein is composed 
of adherends made out of two dissimilar materials, an aluminum sheet 
and a composite sheet, that are bonded to each other with a relatively 
thin, room- temperature curing adhesive. This system is intended to 
model the repair of a cracked, aluminum aircraft component that is 
repaired by bonding a composite sheet to it. Each of the three 
constituents of the system - the metal, the composite, and the adhesive 
exhibits different fatigue behavior and plays an important role in the 
fatigue behavior of the reinforced system. Consequently, to analyze 
the fatigue process in the system, the fatigue behavior of each of the 
constituents needs to be understood. In the following sections the 
fatigue behavior of each of the constituents will be discussed and 
analysis methods formulated. 

Fatigue of Metals 

As pointed out by Erodogan (1968), the fatigue process in metals 
occurs in three different stages: crack initiation, stable crack 

growth, and unstable crack growth (fracture). Current aircraft design 
methods focus on the latter two stages of the fatigue process by using 
a "Damage Tolerant Design Philosophy" (military specification 
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MIL-A-83444). This philosophy, as far as fatigue damage accumulation is 
concerned, admits that initial flaws such as cracks, exist in aircraft 
components that are fatigue critical, i.e. may fail under cyclic loading. 
But, the philosophy also assumes that these initial cracks grow stably 
and can be detected during periodic inspections before they reach a 
critical crack length. Once the damage is detected, it can be repaired. 
Hence, the validity of this philosophy rests on the accurate prediction 
of crack growth rate and critical crack length. Fracture mechanics 
theory can be used to determine both crack growth rate and critical 
crack lengths. 

Fracture mechanics theory was conceived when Griffith (1921) 
related fracture to an energy balance as the crack extended. In 1957 
Irwin related the stress state at the crack tip to fracture. A 
schematic of a crack tip and equations for the stresses very close to 
it (Sih and Liebowitz 1968) are shown on figure 1 (additional terms not 
shown in the equations have a negligible effect on the stress state 
near the crack tip). As may be seen from the equations, as the distance 
from the crack tip, r, approaches zero the stresses become infinite. 
Consequently, at the crack tip where the stresses are infinite a 
singularity exists.^ The coefficients of the stress distributions, 
kj and k 2 , are the stress intensity factors which are used 
extensively in fracture mechanics. The Mode I stress intensity 
designated by kj is associated with the stresses that deform 
the crack surfaces symmetrically with respect to the original plane of 

^In reality infinite stresses cannot exist in the material and 
local yielding of the material occurs. This local yielding is ignored 
in linear elastic fracture mechanics. 



crac 



ki 0 0 30 
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Fig. 1. Stress distribution near crack tip 
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the crack while the Mode II stress intensity designated by kz is 
associated with stresses that cause shear displacements between the 
crack surfaces. In 1960 Sanders showed that the stress intensities 
were related to the strain energy release as the crack extended. Hence, 
the stress state near the crack tip was related to the Griffith theory 
of fracture, and the foundation of fracture mechanics was formed. 

The stress intensities can be determined from complex stress 
functions determined from the theory of elasticity (Sih and Liebowitz 
1968) as 


ki - ika = 2/2 lim {/z - a $ (z)} 

z ^ a 


( 1 ) 


where 


z = X + iy and i = /H” 

and X and y are the cartesian coordinates and $(z) is the complex 
stress function as developed by Muskhelishvili (1975) that satisfies 
the equations (plane strain or stress) 


+ Oy = 2[$(z) + $(z)] 


( 2 ) 


L'ia - a + a, = 2[z<l>'(z) + ')^(z)] 
xy X y LX/ X / j 


(3) 


where a ,o and a are the stresses in the cartesian coordinates 
x’ y xy 

and 'V(z) is another stress function. The two stress functions are 
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functions of both the loading conditions and the configuration of the 

body and will be discussed in detail in later chapters. 

The stress intensities can be related to both crack propagation 

rates and critical crack lengths. On the basis of the Griffith theory 

of fracture, a critical value of the strain energy release can be found 

and hence, according to Sanders (I960) a critical value of the stress 

intensity can be found. For the material used in this study, 2024-T3 

aluminum, Hudson (1969) showed that the critical value for k is 

ic 

56,000 psi-in’’. Hence, with the use of equation (1) and the appropriate 
stress functions, the fracture can be predicted. 

Cyclic crack growth rates were related to the stress intensity 
by Paris (1961) by the empirical formula 

da/dN = C(k )** 

1 

where da/dN is the crack propagation rate, C is am empirical 

constant and k indicates the stress intensity range during cyclic 
1 

loading. Forman et al (1967) improved this equation by including the 
critical stress intensity k^^ and the stress ratio R, which is 
the ratio of the minimum to maximum stress in the loading cycle, 
in the empirical formula 

da/dN = (4) 

(l-R)k,^ - k 

IC 1 

u 

where c and n are empirical constants and k,_ = 56,000 psi-in . 
11 
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For 2024-T3 aluminum Hudson (1969) showed that the constants and 
n^ have the values 

= 3.22 X 10"^^ 
n^ = 3.38 

As the previous discussion implies, once the stress functions for 
a cracked body are known, the stress intensities can be calculated. 

With the stress intensities both the crack propagation rate and critical 
crack length can be predicted. The crack propagation rate and critical 
crack length can be used in a Damage Tolerant Design Philosophy to 
predict life of aircraft components. The life is predicted by first 
assuming that the structure contains cracks. The lengths of these 
cracks are assumed to be the largest crack detected in the structure 
or the largest crack which can be overlooked due to the resolution 
of the inspection technique. Then, by using the assumed or detected 
crack length and fracture mechanics theory, the number of load cycles 
to fracture can be predicted. On the basis of these calculations, 
inspection intervals are determined to assure that cracks can be 
detected and repaired before they reach a critical length. 

Fatigue of Bonded Systems 

To perform a realistic fatigue analysis of the reinforced system, 
the fatigue behavior of the adhesive in situ, herein called the bond, 
must be characterized. Several researchers have shown that the bond 
deteriorates when subjected to a cyclic load. Within this disserta- 
tion this deterioration will be called debonding. Hoffman and June 
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(1973) studied debonding by recording the debond propagation as a func- 
tion of applied load cycles. They showed that a myriad of factors 
such as type of adhesive, adherend and adhesive thickness, method of 
curing, and aging all affect the debonding. Roderick et al (1976) 
showed that debonding could occur as failure within the adhesive as a 
cohesive failure, between the adhesive and an adherend as an adhesive 
failure, or in the composite material. Because of the variety of 
failure modes, the analysis of the debonding is difficult. The most 
progress in analysis of debonding appears to stem from the energy 
approach developed by Griffith (1921). 

The first application of the energy approach appears to be by 
Rippling et al (1964) in the study of fracture toughness of bonded 
joints. Since Rippling's paper, Mostovoy and Ripling have published 
several other papers on fracture toughness of bonds: Mostovoy and 

Rippling 1966, Mostovoy et al 1967, Mostovoy and Rippling 1971. However, 
a correlation between the fracture energy and the stress state near 
the debond tip has not been made in the bonded systems. Wang et al 
(1976) showed that a primary reason for the lack of correlation appears 
to be the development of large regions of plastic yielding in the 
adhesive. Hence, linear elastic fracture mechanics based on small 
yield zones and stress intensities at a crack tip do not appear 
applicable to bonded systems. 

However, by applying an energy approach, Roderick et al (1975) 
showed that the debond propagation rates can be correlated for specimens 
with different thickness adherends with a Paris (1961) type equation 
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db/dN = c^iGf^ (5) 

where both and are empirical constants for a specific bond 
system and G is the strain energy released as the debond extends. 

As shown by Roderick et al (1976), the parameters c^ and 
r \2 vary for different bonded systems. The bond system used in this 
dissertation was 2024-T3 aluminum bonded to unidirectional boron/epoxy 
with a room temperature curing adhesive. Shell EA-934. For this system 
the empirical constants were determined by methods discussed in 
Appendix A and were found as 

c^ = 3.158 X 10"^ 

nj = 3.616 

With these constants, a value for the strain energy release rate, G, 
and equation (5), the debond growth rate can be predicted. The 
calculation of G for debonding in the reinforced system will be 
discussed in detail in Chapter V. 

Fatigue of Composite Materials 

The term "composite" may refer to a myriad of systems composed 
of a wide spectrum of different types of fibers and matrices. Further- 
more, each system may have widely different fatigue characteristics 
depending upon the fiber orientations, stacking sequences, and loading 
conditions. Durchlaub and Freeman (1974) showed that fatigue damage 
in composites may occur perpendicular to, parallel to, or at an 
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angle to the loading axis depending upon the fiber orientation. Foye 
and Baker (1970) showed that the lives of composite laminates could 
vary as much as an order of magnitude by changing their stacking 
sequences. Reifsnider et al (1974) showed that changing the frequency 
of the applied cyclic load affects both the mode of failure and the 
fatigue life As evident from these observations, the fatigue behavior 
of composite materials is complex. 

Currently the understanding of the fatigue process in composites 
appears primitive although some progress in developing an understanding 
has been made. As pointed out by Salkind (1973), fatigue failure in 
composites can occur in different failure modes such as matrix cracking, 
delamination, and fiber fracture. Also, evidence exists that suggests 
that the fatigue process is a result of primarily matrix deterioration 
(Roderick and Whitcomb 1977). If matrix deterioration is the primary 
cause of fatigue in composites, then the various failure modes could 
be explained in terms of different stress states in the matrix depend- 
ing upon the fiber orientation and stacking sequences of a specific 
laminate. Hence, those laminates in which the matrix is highly stressed 
would most likely degrade under cyclic load while those laminates in 
which the matrix is lightly stressed would not. 

Following this line of thought, composites in which the fibers 
transmit the load, fiber controlled composites, would have long fatigue 
lives while those in which the matrix transmits the load, matrix 
controlled composites, would have short lives. An example of a fiber 
controlled composite is a unidirectional laminate loaded along the axis. 
On the other hand, an example of a matrix controlled laminate is one 
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in which the fibers are at 45° angles to the loading axis. As shown by 
Durchlaub and Freeman (1974), the matrix controlled laminate does degrade 
rapidly under cyclic loading while the fiber controlled laminate does not. 

Two basic approaches to predict the diverse fatigue behavior of 

composite materials are currently being developed by researchers. In 

the first approach, the laminate behavior is described by a statistical 

model (Halpin et al 1972) that relates the static strength and fatigue 

life distributions by assuming that the residual strength of the 

laminates degrades monotonical 1y (Yang and Liu 1977). Because this 

method is based on experimental results, it can be applied easily. 

However, it does not apply to laminates whose residual strength does 

2 

not decrease monotonical ly with applied load cycles. Also, this 
method requires extensive testing every time the stacking sequence or 
fiber orientation changes. 

The second approach as developed by McLaughlin et al (1975) 
couples basic fatigue data on the laminae level with a stress analysis 
to predict both the mode of fatigue failure and the fatigue lives of 
laminates. Because this approach is based on laminae data rather than 
laminate data, it can be used to predict the fatigue behavior of lami- 
nates with different stacking sequences and fiber orientations without 
extensive testing. The major drawback to this approach is its 

2 

Durchlaub and Freeman (1974) showed that the residual strength of 
notched laminates could increase after fatigue loading. 

3 

The analysis originally proposed by McLaughlin et al did not 
consider interlaminar stresses and therefore could not account for 
changes in stacking sequences, but incorporation of interlaminar stress- 
es into the analysis has been done and will be shown in a NASA 
contractors report released in 1978. 
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complexity in attempting to develop simple, realistic stress analyses 
and a failure criterion on the laminae level. 

The state of the art of fatigue analysis, in the author's opinion, 
is still in the early stages of development and not yet capable of 
reliable life predictions for general laminates. As a consequence 
the fatigue behavior of the unidirectional boron/epoxy used in the 
present study cannot be described by relatively simple fatigue analyses 
as was the case for the cracked metal sheet and the bond system. Accord- 
ingly, the fatigue behavior will be determined solely by experimental 
data. 

Shockey et al (1970) showed that unidirectional boron/epoxy 
laminates that were loaded along the fiber axis had an average ultimate 
tensile strength of 193 ksi; when these laminates were cycled under 
constant amplitude cyclic loading with R = 0.1, they retained 
73 percent of their ultimate tensile strength after 10^ applied load 
cycles. Consequently, in an attempt to prevent failure in the 
unidirectional boron/epoxy, stress along the fiber axis (based on 
laminate analysis) was kept below 140 ksi. 

Having discussed the fatigue behavior of the constituents of the 
reinforced system in this chapter, the next chapter deals with the 
fatigue behavior of the constituents in situ in the reinforced system. 
Hence, the next chapter discusses fatigue tests of reinforced systems. 



CHAPTER III 


FATIGUE TESTS OF REINFORCED SYSTEMS 

To determine the fatigue behavior of the reinforced system, two 
large panels were manufactured and tested. The panels shown on 
figure 2 were made of 8 x 24 inch sheets of 2024-T3 aluminum and 
unidirectional boron/epoxy. EA-934 room temperature curing adhesive 
was used to join the sheets with the bonding process described in 
Appendix A. The primary difference betv/een the panels labeled A and B 
on figure 2 was the thickness of the metal and composite adherends. To 
simulate a crack, the metal adherend contained a through-the-thickness 
narrow slit 0.01 inch wide and 2 inches long. The slit, which was made 
by an electrical discharge process, was centered along the horizontal 
centerline of the panels. In both panels the fibers of the unidirec- 
tional composite run parallel to the longitudinal axis of the panels. 

The panels were tested in a 300,000 pound load capacity servo- 
hydraulic fatigue machine. Both panels were tested under constant 
amplitude loading with R, the ratio of the minimum to maximum stress 

4 

in the load cycle, equal to 0.01 at a test frequency of 2.5 Hz. 

For the fatigue tests of both panels the distance between test 
machine grips was 16 inches. The maximum loads applied to each panel 

4 

The test frequency was limited to 2.5 Hz instead of the 10 Hz 
used to characterize the debond behavior in Appendix A because of 
test machine limitations. 


16 
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y 



thickness 


Panel 

metal 

composite 


t 

t^ 


m 

c 

A 

0.156 in 

0.028 in 

B 

0.051 in 

0.048 in 


Fig. 2. Panel configuration 
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during the fatigue tests and the corresponding stresses in the adherends 
calculated from membrane laminate theory (see Appendix C) are shown 
below. 


Panel 

Maximum Load 
lbs 

Stress, psi 
metal composite 

A 

37,600 

19,600 

58,600 

B 

22,500 

14,600 

43,100 


During the fatigue tests, crack lengths were measured periodically 
with an optical microscope. Table 1 shows the measured crack lengths 
and applied load cycles for both tests. The crack lengths are plotted 
against the applied load cycles on figure 3. Note that on the figure 
the abscissa is logarithmic and the ordinate starts at the initial half- 
crack length of a = 1.0 inch. The crack propagation rates for the 
panels are the slopes of the curves shov^n on figure 3. These rates 
are tabulated in Table 1 and plotted against the half-crack length on 
figure 4. As evident from figure 4, the crack growth rate is about 
two orders of magnitude larger in Panel A than in Panel B. 

The crack propagation rates in these panels is a function of 
debonding between the adherends. If the adherends were completely 
debonded the crack propagation rate would be much larger than if no 
debonding occurred. To investigate the effect of debond size on 
crack propagation rate, the test panels were examined with an ultrasonic 
C-scan (details of the C-scan method are discussed in detail by 
McMaster 1963) after the half-crack length grew to 1.0 inch. Figure 5 
shows the C-scans of the panels. On the figure the dark parts of the 
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TABLE 1 


CRACK 

LENGTHS 

AND CRACK PROPAGATION 

RATES 




Panel A 

Panel 

B 

half-crack 

cycles 

crack 

cycles 

crack 

length 


propagation 


propagation 



rate 


rate 

a, in 

*N 

da/dN 

*N 

da/dN 



X 10-5 


X 10-5 

1 .05 



44,800 

0.236 

1.10 

1 ,570 


65,940 

0.312 

1.15 


7.57 

81,975 

0.217 

1.20 

3,220 


105,000 

0.33 

1.25 


9.02 

120,130 

0.369 

1.30 

4,440 


133,675 

0.311 

1.35 


10.8 

149,740 

0.328 

1.40 

5,540 


164,980 

0.339 

1.45 


11.6 

179,695 

0.321 

1.50 

6,790 


195,250 

0.277 

1.55 


12.0 

213,310 

0.299 

1.60 

7,790 


230,025 

0.275 

1.65 


12.2 

248,180 


1.70 

8,690 


.... 


1.75 


13.0 

279,795 

0.295 

1.80 

9,690 


296,765 

0.323 

1.85 


13.0 

312,260 

0.339 

1.90 

10,680 


326,990 

0.312 

1.95 



342,990 


2.00 



. . . . 



*N - Instead of listing the number of cycles that caused crack 
growth at both crack tips, the average number of cycles is given. 



half-crack length 
a, inches 



N, applied load cycles 

Fig. 3. Crack growth versus applied load cycles 
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Debonding 

from 

grips 


Specimen B 


Specimen A 


Ultrasonic C-scans of fatigued specimens 
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C-scans are areas where debonding has occurred. As evident from the 
figure debonding occurred over an elliptical area. Ratwani (1977) 
has observed similar elliptical debonds in metal laminates. The major 
axes of the elliptical debonds were nearly equal to the crack length 
in the metal adherend of the panels. As measured from figure 5 the 
minor axes of the debonds were 3.0 inches and 0.50 inch for Panels A 
and B respectively. However, when the C-scans shown on figure 5 were 
made, they were distorted along the longitudinal axes of the panels. By 
taking into account this distortion, the minor axes of Panels A and 
B were found to be 4.0 and 0.67 inches respectively. With these 
corrected values of debond length along the minor axes, the ratios of 
minor to major axes of the debonds, F = b/a, were found as 1.0 and 
0.14 for Panels A and B respectively. 

The experiments discussed in this chapter showed that the fatigue 
in reinforced systems occurs as col linear crack growth and debond 
growth over an approximately elliptical area. In Chapters IV and V 
an analysis will be developed that can model this observed behavior. 

In Chapter VI the analysis will be verified by comparing the experimental 
results of this chapter with results of the analysis. 



CHAPTER IV 


STRESS ANALYSIS 

As shown in the previous chapter, under cyclic loads the rein- 
forced system exhibits both crack and debond growth. Intuitively, 
the rate at which the debond and crack propagates is a function of the 
stress state and level in the system. Consequently, to predict these 
rates a realistic stress analysis is required. For the stress analysis 
to be realistic it must predict stresses in the adhesive as well as 
in the adherends of the system. Because adhesives typically exhibit 
nonlinear behavior (Hughes and Rutherford 1968), the stress analysis 
must include nonlinear behavior of the adhesive. The first step in 
development of a realistic, nonlinear, elastic stress analysis is the 
formulation of a linear elastic stress analysis. 

Formulation of Linear Elastic Solution 

To formulate an elastic solution, the reinforced system shown in 
figure 6 was considered. As shown in figure 6, the system consists of 
a cracked metal sheet bonded to a composite sheet with an elliptical 
debond between the sheets. The system was subjected to a remote 
stress, s. A rigorous stress analysis of this system requires a 
three-dimensional formulation. Although a general, exact solution is 
not available, finite element or finite difference numerical methods 
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can be employed to obtain a tractable solution. However, these solutions 
are not efficient for analyzing reinforced systems in which the crack 
length and debond area are continually changing. An alternate, simple 
analysis can be developed by assuming that the adherends are in plane 
stress while the adhesive is in pure shear. These assumptions were 
first used in an analysis by Volkerson in 1934. 

The validity of these assumptions for analysis of the reinforced 
system shown in figure 6 was investigated in detail in Appendix B. 

As shown in Appendix B with a simple example, the assumption can lead 
to errors as much as 100 percent in the calculated adhesive stresses 
as compared to more rigorous finite element solutions. Evidently, 
significant shearing deformation occurs in the adherends of the 
reinforced system. The presence of the adherend shearing deformation 
violates Volkerson's assumptions, but as shown in Appendix B an 
effective shear modulus, can be determined and used with the 

assumptions to calculate adhesive shearing stresses within a few percent 
of the finite element results. 

Arin and Erodogan (1972) used Volkerson's assumptions with complex 
variable elasticity theory developed by Muskhelishvili (1953) and 
Lekhnitski (1956) to analyze a system similar to that shown in figure 6. 
The linear elastic stress analysis developed herein basically follows 
the concepts used by Arin and Erodogan, but differs in the formulation 
of the Green's functions used in the elasticity solution, the method of 
numerical integration of the Green's functions, and the domain of 
integration. To develop the stress analysis, the reinforced system 
is free bodied as shown in figure 7 (adhesive layer not shown) 




nforced system 
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using Volkerson's assumptions. In the figure the remote stress, s, 
refers to the applied load over the total cross sectional area of the 
reinforced system. By the use of laminate theory as described in 
Appendix C, the inplane stresses a , a , a , and a (where 

Jr II* J' ^ /\lll 

the first subscript refers to stress direction and the second subscript 
refers to the metal or composite adherend) can be easily calculated. On 
figure 7, f^^ and f^^ indicate shear stresses in the adhesive 

layer. Throughout this dissertation the adhesive shear stresses which 
will be assumed to act as body forces on the adherends will be frequently 
called interlaminar stresses. To form a governing equation, these 
interlaminar stresses were related to the displacement of the adherends 
in the following manner. 

First, the shear strain in the adhesive layer was related to the 
displacement in the adherends by the relations 


Y 


xz 


m 



Y 


yz 



( 6 ) 


where t^^ is the thickness of the adhesive, u and v refer to 
displacements in the x and y direction respectively, and the 
subscripts m and c refer to the metal and composite adherends 
respectively. Next, Hook's law was used to relate the shearing strain 
in the adhesive to the interlaminar stresses as 


T = GY 


(7) 
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Then, by the use of the effective shear modulus and equations (6) and 
(7) the interlaminar stresses can be related to the adherend dis- 
placements as 


’eff 


"eff 


xz 


‘“m - “c’- 


"ad 


zy 


^'m - "c> 


(8) 


"ad 


Equation (8) can be written for several different points in the 
reinforced system to form a system of simultaneous equations. These 
simultaneous equations are developed in detail in subsequent sections 
of this chapter. 

The displacements, u and v, in the adherends were related 
to the inplane adherend stresses and the interlaminar stresses by 
several functions - Fg as 


u 

V 


m 

m 


’'.("xx 

''^(“"xx 


.Olllyy) ^S^^XZ*^yZ^ 

’™yy> * ^<^z-fyz> 


(9) 


u = F (ac ,ac 
c 3 ' xx’ yy 




The displacements. 


> * F,(fxz'fyz 

^ ^e^^xz’^yz 
Fj and Fj, 


) 

) 

in the metal adherend due 


to the remote stress were calculated in terms of two stress functions 
(|)(z) and (o(z) (Muskhelishvili 1975) by the equation 
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2G(Ujj^ + = n<l){z) - u(z) - (z - z)$(z) (10) 


3 - V 

where n = for plane stress, i = z = x + iy, the bar 

1 + V 

over a function or variable denotes the complex conjugate, v 
is Poisson's ratio and 



( 13 ) 


In equations (11) through (13) "a" denotes half the crack length 
in the metal sheet. With the use of equation (10), the functions F 
and F were found to be 
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F (am ,ani ) = Real 
1 XX’ yy' 


ncl)(z) - w(z) - (z - z)$(z) 
2G 


(14) 


F^(am^^,aitiyy) = Img 


in(t)(z) - u(z) - (z - z)$(z) 
2G 


The displacements, Fj and F^, in the composite sheet due 
to the remote stress were found as follows. First, the constitutive 
equations for an orthotropic material were used to relate strainsto 
stresses (Lekhnitskii 1968) as 


ac V 
XX yx 


V 


xy 


acyy. ey=- 


— ac + — 
E XX E 
X y 


ac 


yy 


(15) 


where E and E are the moduli of elasticity in the x and y 
A y 

directions respectively, and v and v are Poisson's ratios. 

xy yx 

Then with the definition of strains as e = and e = It 

X dX j oy 

equations (15) were integrated to find displacements as 


ac V 
XX yx 

u = F,(cJC_,ac„J ={ ac,„> x + hl(y) 


3 ' XX ’ yy ' 


V = F (ac ,ac ) = < 
XX ’ yy' ' 




V 


yy 


ac. 


xy -'yy, 

ac + >y + h2(x) 

,x E 


( 16 ) 
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where hl(y) and h2(x) are arbitrary functions which were set to 
zero because of symmetry considerations. 

The displacements, and Fg, in the metal adherend 
caused by the interlaminar stresses were calculated by assuming 
that the interlaminar stresses acted as body forces on the adherends. 
With this assumption, the displacements were calculated using Green's 
functions in surface integrals as 


F, (f ,f ) = aa + aa 
1' xz’ yz' XX xy 


F,(f ,f ) = aa + aa 
2' xz’ yz^ yx yy 


where 


( 17 ) 


“xx = "s =‘>xx<^’^o>'-<'xz(^o>]''’'o<‘l'o 

( 18 ) 

“yx = "s 8Dyx(z-2o»-fxz<^o>'“''o'^l'o 

“yy ' «'>yy(Z'^o>!-''y2(^o)5<'Vl'o 

and Zq is the location of a point load (see Appendix D) and GD^^, 

GD , GD , and GD are the Green's functions which were 
xy yx yy 

discussed and derived in Appendix D and found as 
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GD/x = CQ{Real(g(z,ZQ) + g(z,ZQ))} 
GDxy = CQ{Real(i{g(z,z^) - g(z,z^j)})} 

GDyx = CQ{Img(g(z,z^) + g(z,z^))} 

GDyy = CQ{Img(i{g(z,ZQ) - g(zJjj))} 


(19) 


where 


g(z,ZQ) = n[Real{XA(z,ZQ)} - XC(z,z^j)] 


+ 0.5[n^XB(z,ZQ) + XB(z,Zq)] + XC(z,Zq) 


+ 2 


“o - ^ 0 ^ 


z - z. 


+ (z - z)B(z.Zjj) 


with 


XA(z,Zq) 


= - Log <- 


zzo - a^ - Kz^)I(z): 
zZq + a^ - I(Zq)I(z)I 


XB(z,Zq) 


= Log 


zZg - a^ - KZjj)I(z) 

r-_^oi 

zz„ + - I(z )I(z) 

0 ' 0' ' ' 

" 1 

N 

1 

O 

N 

1 
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and 


•4z, 


2nz, 


z - z. 


B(z,Zj^) = 2 




fz^ - Z, z^ - z^^ (z - z,)^ (z + z^) 


J 


z„ - z^ \ a - zz„ a + zz„ 
0 0 1 0 0 


2I(z^)I(z) [iz-z^V (z + Zo)'* 


I(z) 


KZq) 


KZo) 


- n- 


Z„2 - z2 72 - z2 
0 0 


I(z) = \/z2 - a2 


^0 


4Gt^Tr(l + n) 


The domain of integration used in equation (18) will be discussed 
in a later section. 

The displacements, F 7 (fy,»^w,) and F (fy,.f,,,)> in the composite 
adherend caused by the interlaminar stresses were found with an 
approach similar to that used for the metal adherend. The functions 


were written as 
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F_{f ,f ) = bb + bb 
7' xz’ yz' vx i 


XX xy 


F (f ,f ) = bb + bb 

S\ XZ’ yZ^ >JyX ^ yy 


where 


( 20 ) 





0 


% “ ”'’xy<^'"k>fy2(\)'‘V>'o 

(21) 

‘>'=yx = ™yx<^k-'*k>fxz<"k><‘>‘o«l'o 

bb^y = ff HDyy(Z|,.w^)fy^(W|,)dx„dy^ 

k = 1,2 

and HD , HD , HD , and HD are the Green's functionsfor the 
XX xy yx yy 

composite adherend derived in Appendix D and found to be 



2Real{pjC11g^ (Zj,Wj) + P2C21g^ (z^,\/2)} 


HD^y = 2Rea1{i[pjCl2g2(Zj,Wj) + P2C22g2 (z^ .w J]} 
HDy^ = 2Real{qiC11gi(Zj,Wj) + q2C21gi(z^,v/^)} 


( 22 ) 


HDyy = 2Real[i[qjC21g2(Zi.Wi) + q2C22g2(z2,w ,)]] 
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Zj = X + SjY, = X + s^y, = Xq + 

w = x„ + s,y^ 

2 0 

where and are complex numbers which are not complex conjugates 
of each other and are roots of the equation (Lekhn tski 1968) 
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Substituting equations (14), (16), (17), and (20) into equation (9) 
and that result into equation (8) yields the governing equation that 
was used to formulate a system of simultaneous equations. 

Numerical Solution 

To solve for the unknown stresses f and f , the domain of 

Ai y i 

integration in equations (18) and (21) was separated into three 
regions as shown on figure 8. Region A on the figure represents 
the portion of the domain where debonding has occurred. In this region 
the interlaminar stresses are zero. Region B on the figure represents 
the portion of the domain where significant interlaminar stresses 
exist. As shown on the figure, this region is divided into smaller 
elements. Region C on the figure represents a portion of the domain 
where the interlaminar stresses are small and can be neglected. The 
size of each of these regions will be discussed further in Chapter VI 
where convergence of the system is investigated. 

The next step in the formulation of the simultaneous equations 
was to assume that the interlaminar stresses were constant over each 
element of region B. With this assumption, the displacements caused 
by the interlaminar stresses, equations (17) and (20), were written 
as 
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y 



Fig. 8. Domain of integration 
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n 

F^(i) = E AA^^(i,j)[-fj^^(j)] + AA^y(i.j)[-fy2(j)] 
J=1 


^ AA^^(i.j)[-f^2(j)] + AAyy(i,j)[-fy2(j)3 

J=1 


F^(i) = E BB^y(i,j)[-f^^(j)] + BB^y(i,j)[fy^(j)] 

J=1 


( 23 ) 


F^(i) = E BBy^(i,j)[fj^^(j)] + BByy(i,j)[fy^(j)] 
J=1 


where n is the number of elements in region B and 


AA^^O.j) = /; GD^^(z,,z„)dx„dy„ 
AA^y(l.j) = ff GDj,^(z,.z„)dXj,dy„ 
AAj,^(7,j) = ff GDy^(z,.z„)dx„dy^ 
AAyyd.j) = ff 

BB^,(i.j) = ff HD^^(z^,.W|^)dx„dy„ 
= ff HD^^(z^,,W|^)dx„dy„ 
= ff HDy,(z^,.w^)dx„dy„ 


BB^^(i.J) = fl HD^j,(z^,,W|^)dx„dyj, 


k = 1,2 
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where the i subscript indicates the point in the x-y plane, z, 
where the governing equations (equations (8)) were evaluated. The 
j subscript when used in the coefficient of the interlaminar stresses, 
f and f , indicates the location of the centroid of the element 

JT 

over which the interlaminar stresses act while the j subscript 
used in the interlaminar stresses indicates the value of these stresses 
acting on element j. The integrals in equations (24) were evaluated 
numerically; the method of integration will be discussed in Chapter VI. 

Substituting equation (14) and (23) in equations (9) and (8) and 
evaluating the result at the centroid of each element of region 
B lead to a system of 2n simultaneous linear equations where n 
is the number of elements in region B as 

AAxx(i.j) + BB^^(i,j) AA^y(i,j) + BB^^(i,j) 

AAy^(i,j) + BBy^(i,j) AAyy(i,j) + BByy(i,j) 


^ad 

+ 

1 

0 


F (1) - F (i) 

1 3 

^eff 

0 

1 




Using Gaussian elimination, this system of linear simultaneous 
equations was solved and yielded values of the unknown interlaminar 
stresses and 
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Nonlinear Solution 

As shown in Appendix B, the adhesive used in the reinforced 
system can exhibit nonlinear stress strain behavior. As shown in the 
appendix, the tensile stress strain curve for the adhesive can be 
approximated by a bilinear stress strain curve with a change of slope 
occurring at about 4200 psi. Because the adhesive in the reinforced 
system is assumed to be in a state of pure shear, the data from the 
tensile stress strain curve must be related to the adhesive in a 
pure shear state. To develop this relationship, a yielding criterion 
is required. 

For simplicity the criterion developed by Von Mises and given 
by Hill (1951) as 


(a - a )^ + (a - a )^ + (a - a 

' W \t\t f ' V/\/ ^ ~9 f ' -7-7 VV' 


XX yy' 


yy zz' 


ZZ XX' 


(26) 




where a is the stress at yielding and is a constant, was used to 
estimate when the adhesive in the reinforced system would exhibit 

c 

nonlinear behavior. For pure tension, as was the case in the bulk 
property test described in Appendix B, equation (26) reduces to 

^Several yield criteria, of which Von Mises' is one of the most 
popular, are available in the literature. 
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•N 



( 27 ) 


For the adhesive in the reinforced system the only stresses present, 

according to the Volkerson assumptions, are a and a , Hence, in 

yz xz 

this case equation (26) reduces to 


(a‘ 


yz 


^ 1 = k ^ 

xz' ^0 


(28) 


Combining equations (27) and (28) by eliminating k^ yields a relation 
between tensile and shear yielding as 



+ 


xz 



/3 


(29) 


For the bulk property tensile tests was found approximately 

at 4200 psi. With this value in equation (29) and the notation for 
interlaminar stresses in the reinforced system an inequality was 
developed as 


f"yz + > 5.88 X 10^ (30) 

Equation (30) was used to estimate the initiation of nonlinear behavior 
of the adhesive in the reinforced system by using f and f 

Jr Z XZ 

from the solution of equation (25). As will be shown in Chapter VII, 
equation (30) predicts that nonlinear behavior of the adhesive will 
occur in many reinforced systems. The adhesive nonlinearity manifests 
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itself in equation (25) as changes in with the magnitude of 

applied load. As shown in Appendix B, the for the adhesive 

can be either 65,000 psi or 36,000 psi. The value G^^^ takes in the 
reinforced system can be determined from equation (30). If equation (30) 
is true then G^^^ = 36,000 psi while if equation (30) is false 
^eff ” psi • 

By the use of inequality (30) and equation (25), the applied 

stress at which G^^^ of the adhesive changes in the reinforced 

system was predicted with the following approach. As shown by equations 

(14) and (16) the right hand side of equation (25) is a function of 

the applied stress, s (om , am , ac , ac are linear functions 

XX yy XX yy 

of s). Hence, the solution of equation (25) was written in terms 
of the unit solution and the applied stress, s, as 

" ^^xz^'^^unit’ ” ^^yz^"^^unit (31) 

where s is the remote applied stress and f^z^^^unit ^yz^^^unit 

are the solutions of equation (25) with an applied stress of s = 1. 
Substituting equations (31) into equation (30) and solving for s 
for each of the elements of region B, lead to values of the remote 
stress s(j) which cause a change in G^^^ for element j as 
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The smallest value of s(j) for n elements of region B is the value 
of the applied stress, s^^^, which causes a change in in 

element k. At this value of s^^^ all elements except element k 
have = 65,000, while element k has = 36,000. 

If the value of s^*^ is greater than the maximum remote 
applied stress, then the solution is completely elastic and 

from equation (31) the interlaminar stresses in the system are given 
as 


^max^xz^"^ ^unit 

(33) 

fyz(j) " ^max^yz^'^^unit 

However, if s^^^ is less than s„.„ then yielding has occurred.^ 

iTiaX 

At the yield point the stress in all of the elements is given by 

(34) 

9yz‘*(j) ' 

where the unit(i) indicates that the unit stresses were obtained from 
an elastic solution where all had the same shear modulus of G^^^ = 

^Yielding refers to a change in G^^^ in the adhesive of the 
reinforced system. 
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65,000 psi. Once yielding has occurred in element k the shear 
modulus of element k takes on the secondary value = 36,000. 

Consequently, for the next increment of applied stress the 
governing equations (25) must be modified as 


1 


[A] 


+ 




eff 


G2 


eff 


ff 


" 

xz 

= 

*^unit 

f 



yz 





■ 


(35) 


where 


[A] = 


AA^^d.j) + 


AAy^(i,j) + BBy^(i,j) 


AA^y(i,j) + BB^y(i,j) 
AAyy(i,j) + BByy(i,j) 


Equation (35) was then used to find a new unit solution after 
element k had yielded. The new unit values, ^unit(z) 

^yz^^ ^unit(2) ’ equation (31) and added to equations 

(34) to give the stress in each element after the second load increment 
as 

f„U) = gJ'’(j) + 


(36) 
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/ 2 ^ 

Substituting equations (36) into equation (30) and solving for s' ' 

{ 2 \ 

yields the value of s' ' that causesthe next element to yield 
as 


(^) 


-B^ + V 
0 Y 0 


- 


2A, 


(37) 


where 


'o “ |^xz^'^^unit(2^ |*'yz^'^^unit(2)| 

' 4gJ‘>(j)f„(j)„„,t(0 * 9yz‘’(J>V‘^’uni 


unit( 2 ) 


( 2 \ 

If the value of s' ' > s then the stresses after the section 

max 

increment of load are given by 


' <=max ‘ fxz'j>unU( = ) ^ 9xz‘’<J> 


- s''*) t g„i‘*(j) 


^yz 


yz''^ 'unit( 2 ) =yz 


( 38 ) 
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( 2 ) 

and the total number of load increments is 2. But if s' ' is less 

than s„,„ then the stresses after the second load increment are 
max 

given by 


(38) 


The entire process is repeated until the sum of all the load increments 

equals The final value of the interlaminar stresses are the final 

max 

values of and gy^(j) given as 


fx,(j) = gx'”>(0) 


(40) 


where m is the total number of load increments. 

To obtain the final values for the interlaminar stresses, equation 

(40), the system of simultaneous equations (equations (35)) must be 

solved for each load increment. To minimize computational effort, a 

Gauss Siedel method discussed by McCracken and Dorn (1968) was used to 

solve equation (35) after the first unit vector was found by Gaussian 

elimination. By the use of unit vectors, ^unit(k) 

^yz^"^ ^unit(k) ’ ^ iteration as initial estimates for the k + 1 

unit vectors, the Gauss Siedel method rapidly converges. Consequently, 

the method is very efficient in solving for the unit stresses in 
successive load increments. 
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In this chapter an analysis was developed to determine inter- 
laminar stresses in the reinforced system. In the next chapter these 
interlaminar stresses will be used to determine crack and debond 
growth rates in the reinforced system. The accuracy of the analysis 
will be discussed in Chapter VI. 



CHAPTER V 


FATIGUE ANALYSIS 

As shown in Chapter II, fatigue damage in the reinforced system 
occurred as crack and debond growth. Accordingly, an adequate fatigue 
analysis should predict both types of growth. To the author's 
knowledge, the only analysis to date that attempts to model this behavior 
was developed by Ratwani (1977). Ratwani (1977) analyzes crack and 
debond growth in a cracked metal sheet reinforced with an uncracked 
metal sheet by using the elastic analysis developed by Erdogan and 
Arin (1972) and a maximum strain criterion to predict debond growth. 

In contrast, herein, the nonlinear elastic stress analysis developed 
in the previous chapter was coupled with a debond propagation equation 
based on calculated strain energy release rates. 

Crack Growth Rate 

To use the crack propagation equation (equation (4)), the stress 
intensity must be determined. The stress intensity can be related to 
the two stress distributions discussed in Chapter IV: the remote 

inplane stresses in the reinforced system and the interlaminar shear 
stresses which act near the debond front. The stress intensity can be 
found by superimposing the stress intensities from the two stress 
distributions. 
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The stress intensity produced by the remote stresses can be 
calculated by substituting the stress function, $(z), for this 
loading case (equation (13)) into equation (1) which related the 
stress function to the stress intensity as 



which results in 


k = 0 


(42) 


The stress intensity produced by the interlaminar shearing 
stresses, fy,(j) and f (j), was found in the following manner. 
First, the stress intensity for four point loads acting on a cracked 
sheet was found by substituting the stress function 


<I>(z,Zq) = SB(z,Zq) + SB(z,Zq) 


(D.43) 


derived in Appendix D into equation (1) to yield 
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kj - ikj = 2/2 litti ^ - a 

z -> a 


SB(z,Zq) + SB(z,Zq) 


(43) 


Taking the limit of equation (43) and combining coefficients of the 
X and Y load components results in 


kj - ik^ = X[XK(Zq) + XK(Zp)] 

+ iY[XK(z^^) - XK(z^)] 


which leads to 


k^ = 2Real[XK(z^)]X - 2Img[XK(zQ)]Y 

kj = 0 


(44) 


where 


/a 


XK(2„)= 


^0^0 - “o' * 




nl(z„) 




(45) 


Then, with the use of the coefficient of the X and Y load 
components as Green's functions for the stress intensities, the stress 



intensity caused by the interlaminar stress, and f (j), 

Afc Jr 

can be found as 


n 

ki = 2 Z XK{Zq) dx^dy^ 

i=l 

(46) 

+ fy2^J^-^^"’9(XK(zQ)dxQdy^] 


The domain of integration for equation (46) is, of course, the same as 
was discussed in the previous chapter. The method of integration will 
be discussed in the next chapter. 

Debond Growth Rate 

To use the debond propagation equation (equation (5)), the 
strain energy release rate must be determined for points along the 
debond front. A rigorous determination of the strain energy release 
rate is difficult and beyond the scope of this effort. However, an 
approximate solution is developed in the following paragraphs. 

A freebody of a strip was taken from the longitudinal centerline 
of the reinforced system as shown in figure 9. The strain energy 
release rate for the strip was approximately calculated with the J 
integral developed by Rice and given in Liebowitz (1968) as 
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where a., is the stress, e. is the strain, T is the surface 

I J I J I 

traction, u^. is the displacement and r is any contour that 
contains the debond front and does not pass through a plastic region 
of material. Equation (47) can be written in terms of cartesian 
coordinates as (Yoder and Griffis 1974) 


9v 


9wT 




yy 


- a. 


9v 


9y 

9w 


yz 


> dz 


3y 


(48) 


+ <a 


- + a — > dy 

^ 9y 9y 


where is the strain energy density and v and w are displacements 

in y and z directions respectively. 

To apply integral (47) to the freebody shown in figure 9, a path 

of integration shown in figure 10 was used to analyze the energy 

release rate. The path surrounds the debond front on which no 

stresses or traction act, follows the bond line in the metal adherend 

on which the interlaminar stresses act, crosses the adhesive which is 

assumed to have insigificant stresses, and follows the bond line in the 

composite adherend on which the interlaminar stresses act. The 

9v„ 9v 

n c 

— is the strain in the metal adherend, and the — is the 

9y 9y 


Strain in the composite adherend. With the use of this contour, the 
value of G from equation (48) can be written in terms of the inter- 
laminar stresses, f^^ and f^^, and the strain in the adherends, 

V 
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Fig. 10. Contour for determination of G for strip 
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With the use of discrete values of the interlaminar stresses and strains 
in the adherends, the strain energy release rate G can be approximated 
as 

P 

1=1 

where represents the length of the discrete elements, i the index 
for the elements in the strip, and p the number of elements in the 
strip. 

The interlaminar stresses, fy^(k), were determined from 
equation (40). However, the strains in the adherends still need to 
be determined. The strains in the adherends were determined from the 
stresses in the adherends according to laminate theory as shown in 
Appendix C as 



metal composite 


( 51 ) 
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where matrices C and D are given by equations (C.4) and (C.5) 
of Appendix C. 

Before equation (51) was used to calculate the strain in the 
adherends, the stresses in the adherends were determined. 

As shown by Muskelishvili (1975) the stresses in the metal 
adherend can be expressed in terms of two functions, $(z) and Sl(z), 
as 

= Real{3$(z) - Q(z) - (z - z)o' (z)} (52) 

Oy = Real{$(z) + ^{z) + (z - z)$' (z)} (53) 

a = Imag{$(z) + a(z) + (z - z)$' (z)} (54) 

Ay 

Equations (52) - (54) were used to determine the stresses in the metal 
adherends caused by both inplane remote stresses and interlaminar 
stresses. 

For remote inplane stresses, 4>(z) is given by equation (13) 
as 

^(z 

<I>'(z) was found by taking the derivative of $(z) as 

^yy ^ 

2 ( z ^ - a ^)'/2 



$i(z) = - 


(55) 
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and f2(z) was found by differentiating equation (12) as 



For the interlaminar stresses, $(z) and 'i'(z) were determined 
and substituted in equations (52) - (54) (see Appendix D). The result, 
stress in the metal caused by interlaminar stresses, was found as 

n 

j=' 

n 

= -E + GSyj,fy^(j))dx„dy^ (58) 

j=l 

n 

<^xy = + GS(^y)yfy^(j)}dx^dy^ (59) 

j=l 

where GS^^, GS^^, GS^^^, GS^^. and GS^^yjy are the Green's 

functions for stresses given in Appendix D by equations (D.68) 
through (D.73) respectively. The domain of integration in equations 
(57) - (59) is the area of an element shown in region B of figure 8 
where j denotes the particular element. 

Adding equations (52) - (5?) and (57) - (59) yielded the stress 
in the metal adherend of the reinforced system. This result was then 
used in equation (51) to calculate the strain in the metal adherend. 

f 
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The inplane stresses in the orthotropic adherend caused by the 
remote stresses were calculated with simple laminate theory as shown 
in Appendix C. 

The inplane stresses in the orthotropic adherend caused by the 
interlaminar stresses were found in Appendix F as 

n 

n 

j = l 
n 

”xy= ♦HS,^)yfj,^(j))dx„dy„ (62) 

j=l 

»t.ere HS^,. HS^, HS^^. HS^^. and are tha Green's 

functions for an orthotropic solid sheet given in Appendix E by 
equations (E.23) - (E.28). 

Adding the inplane stresses in the orthotropic adherend caused 
by the remote stresses and the interlaminar stresses and substituting 
the result into equation (51) yielded the strains in the orthotropic 
adherend. 

With the strains in the adherends, equation (50) was used to 
calculate the strain energy release rate for the debond along the 
longitudinal axis of the reinforced system. 
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As indicated in Chapter III, the shape of the debond throughout 
the cyclic tests can be approximated by an ellipse. The equation 
for a general ellipse is given as 

y = b[l - (ijn'*** (63) 

where a is half the crack length and b is the debond length 
measured along the y-axis from the center of the crack. With the use 
of equation (63) the overall debond shape was predicted by calculating 
the half crack length, a, from equation (4) using the stress 
intensity calculated from equation (45) and calculating the debond 
length, b, from equation (5) using the strain energy release rate 
calculated from equation (50). 

Prediction of Crack and Debond Growth 

The analysis developed in this dissertation was programmed on 
a CDC 6600 computer at NASA-LRC. A discussion of the program, a 
sample analysis, and a program listing are given in Appendix F. 

A flow chart of key elements of the program is shown in figure 11. 

In the next chapter the accuracy of the analysis is investigated 
by comparing results of the calculations from the analysis with 
experimental results generated in Chapter III. 
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Fig. 11. Flow chart for fatigue analysis 













CHAPTER VI 


ASSESSMENT OF ACCURACY 

Before the analysis developed in the previous chapters was used 
to study the fatigue behavior of the reinforced systems, the accuracy 
of the analysis was assessed. This assessment was based on analytical 
convergence studies and studies comparing analytical results with 
experimental results obtained in Chapter III. 

Numerical Integration 

A key item in the analysis was the method of integration of the 
Green's functions, In theory the Green's functions for both stress 
and displacement require integration over an infinite domain. Because 
the functions are complicated, a closed form integration is difficult 
if not impossible to perform. Consequently, a numerical solution was 
employed. Two key items in this numerical integration were the domain 
of integration and the method of numerical integration. 

As shown on figure 8, the infinite domain of integration was 
divided into three regions: A, B, and C, Only region B has significant 

interlaminar shear stresses which need be integrated, To perform the 
integration region B was divided into elements as shown on figure 8. 

Each of these elements was bounded by the curves 
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X = Xj X = Xj 

yi = fl(x) Y 2 = fz(x) 


(64) 


where the functions fj and fj are the form of equation (63) and 
reflect the debond shape. For example, the bounding curves for 
elements along the debond front are given as 


y^(x) = b yl - 






yj(x) = (b + db) 


1 ' 




a + da 


(65) 


where a and b are the half crack length and the debond height 
respectively and da and db represent fractions of a and b. 

The numerical method of integration used for each element of 
the region B was a two-dimensional Simpson's integration. Simpson's 
integration was used so that several of the Green's functions given by 
equations (19), (22), (D.68) - (D.73), and (E.23) - (E-28) could be 
integrated simultaniously by using common values of the complex 
functions. In each element the interlaminar stresses were assumed 
to be constant and the Green's functions were integrated by using 9 or 
18 integration points. Nine integration points were used when the 
domain of integration did not contain a singularity while 18 points 
were used when a singularity existed within the element. 
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In the latter case the domain of the element was divided into 
three regions; two of which were analytic and the third, which contained 
the singularity, was nonanalytic. The two analytic regions were 
integrated with the nine point Simpson's integration scheme. The 
nonanalytic region was integrated by separating the Green's functions 
into products of analytic and nonanalytic functions. The analytical 
portions were expanded in a Taylor series and only the first terms, 
which were constants, were retained. The singular portion was integrated 
analytically in the principal value sense (Hilderbrand 1950). The 
product of the first term of the Taylor series and the principal 
value resulted in an approximate integral value in the nonanalytic 
region. In general the value of the integral in the nonanalytic region 
of the element was small in comparison with values in the two analytic 
regions. 

The size of region B shown in figure 8 was determined iteratively 
by starting with a small region and increasing its size until no 
changes occurred in the interlaminar stresses, stress intensity or 
strain energy release rate. As an example. Panel B (under a load of 
22,500 pounds) discussed in Chapter III was modeled as shown in figure 
12 with an initial crack length of one inch and no debond (b = 0). As 
shown in the figure, region B was assumed to be nearly rectangular 

with n elements along the x-axis and n elements along the y-axis. 

X y 

The elements have a width of 0.25 inch and a height of 0.20 inch. The 
analysis was performed for values of n^ ranging from four to six and 
values of n^ ranging from one to six. Test conditions for each of 
these possible combinations were notated as 
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Fig. 12. Domain B for convergence analysis 
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5 
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conditions 

1 

1 

7 

13 

2 

2 

8 

14 

3 

3 

9 

15 

4 

4 


16 

5 

5 

11 

17 

6 

6 

12 

18 


Figure 13 shows the values of interlaminar stresses along the 

longitudinal centerline for the different conditions. As evident 

from the figure, values of the interlaminar stresses have converged 

for values of n = 4 and n = 3. Also both the stress intensities 
X y 

and the strain energy release rates have converged for these values. 
Thus, for this sample case the domain of region B is about 1 x 0.6 inch. 
Similar analyses showed that the length of the domain of region B 
along the x-axis is typically the length of the crack. However, as the 
adherend thickness or shear modulus of the adhesive changes, the 
extent of region B in the y-direction also changes. 

With the use of the one-dimensional analysis discussed in 
Appendix B, the extent of region B in the y-direction was estimated for 
different adherend thicknesses and adhesive moduli in the following 
manner. A strip was taken from along the longitudinal centerline of the 











interlaminar shear stresses, psi 
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Fig. 13. Convergence of interlaminar shear stresses 
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reinforced system as shown in figure 9. Because of the crack the 
load in this strip must be transferred either to the adjacent metal 
via inplane shear stresses or to the composite via interlaminar stresses. 
If all the load were transferred via the interlaminar stresses, the 
interlaminar stresses could be higher and extend over a greater area 
than if the inplane stresses were considered. Consequently, the boundary 
of the interlaminar stress region calculated from the one-dimensional 
analysis, which only considers interlaminar stresses, would be an 
upper bound. 

From the one-dimensional analysis the shearing stresses are 
found by equation (B-12) as 

^(y) = (B.12) 

where 

1 ] ^"^ad 

— + } K = 

Em VE.1 t;,HE9*^ 

m c cj ad 2 

Examination of equation (B.12) revealed that the T(y) is a maximum at 
y = 0. The shear stresses were assumed to be negligible when they 
were smaller than 5 percent of the maximum value. In equation form, 
the relationship between y and 95 percent of the maximum inter- 
laminar shear stress was expressed as 
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0.05 K = K 
0 0 

Solving for y led to 

ln(0.05) 

y = z— (66) 

/a 

the distance at which the shear stresses are 95 percent of their 
maximum value. As evident from equation (66) the y distance is a 
function of the thicknesses and material properties of the adherends 
and the adhesive. For the adhesive used in the reinforced system of 
Panels A and B two effective shear moduli were found as 65,000 psi 
prior to yielding and 26,100 psi after yielding (see Appendix B). 

With the use of the latter of these two values to give a conservative 
bounds, the y distance for region B of Panel B (see figure 2) as 
calculated by equation (66) was found as y = 0.60 in. The estimated 
value agrees well with convergence study results shown on figure 13. 

The same logic applies to both bonded (b = 0) and debonded systems 
(b > 0). Consequently, the domain of region B used in the integration 
of the Green's functions was determined by the crack length and 
equation (66). 

Once the domain of region B was determined, the effect of mesh 
refinement within region B was investigated. Figure 14 shows a 
comparison of interlaminar stresses calculated (again for Panel B) 
using two different mesh sizes in region B. As evident from the figure. 




distance along the y axis, inch 


Fig. 14. The effect of mesh si 



on interlaminar stresses 
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the interlaminar stresses calculated with the fine mesh can deviate 
as much as 16 percent from stresses calculated with the coarse mesh. The 
stress intensities and strain energy release rates calculated using 
the two different meshes were found as 

k 6 

1 

coarse mesh 0.281 1.29 

fine mesh 0.280 1 .36 

The values of and G are not nearly as sensitive to changes in mesh 
size as the interlaminar stresses were. The reason for the insensitivity 
is most likely because both k^ and G are obtained by integration 
schemes that smooth out the effect of local approximations in the 
interlaminar stresses. Consequently, the grid can be rather coarse 
and still accurately predict both k^ and G. In contrast the values 
of the interlaminar stresses require a finer mesh for accurate values. 
Because k^ and G are used to predict the fatigue behavior of the 
reinforced system, a relatively coarse mesh was used in the analysis 
without loss of accuracy. 

Accuracy of the Analysis 

To ascertain the accuracy of the analysis, calculated values of 
stresses in the adherends, the stress intensities, and the crack 
propagation rates and debond sizes were compared to experimental 
results on Panels A and B shown in Chapter III. To compare calculated 
and experimental values of stresses in the adherends and stress 
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intensities. Panels A and B were modeled when the half crack length, a, 
reached 2 inches as shown on figure 15. The debond sizes were obtained 
from the C-scans of the specimens shown on figure 5. The height of the 
grid from the edge of the debond was determined from equation (66). 

The crosshatched elements on the figure indicate elements in which the 
adhesive has changed modulus (yielded) during each applied load cycle. 
The effect of the modulus change on and G will be discussed in 
the next chapter. 

The values of the crack growth rate for Panels A and B were 
calculated for a half-crack length of 2 inches, the debond sizes 
observed in figure 5, and the applied loads shown on page 18. 

The calculated values and the experimental values, from Chapter II, 
of the crack propagation rates are 



da/dN, 

in/cycle 

Panel 

calculated 

experimental 

A 

2.85 X 10"^ 

1.30 X 10’^ 

B 

2.26 X 10"^ 

3.12 X 10“^ 


The difference between the calculated and experimental rates are 
within the scatter of the predicted rates for unreinforced metal 
sheets (Figge and Newman 1967). Hence, the analysis accurately 
predicts the crack growth rate if the crack length and debond are 
known. 

As a further check, the analysis was used to predict the crack 
growth in Panels A and B. Figure 16 shows the calculated and 




Panel A Panel B 


CO 


Fig. 15. Mesh for analysis of Panels A and B 


half-crack length, inches 



applied load cycles 


Fig. 16. Experimental and calculated crack lengths as a function of applied load cycles 
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experimental half-crack lengths plotted against the number of applied 
load cycles. For both panels the analysis gave a conservative prediction 
of the crack growth. The calculated and experimental half-crack lengths 
agree within a factor of 2. This deviation is within the scatter of 
crack length prediction of unreinforced metals. Hence, the analysis 
appears to accurately predict the crack length as a function of 
applied load cycles. 

On figure 17 the calculated and experimental crack propagation 
rates were plotted against the half-crack length. The calculated 
crack propagation rates were within a factor of 2 of the experimental 
rates. The largest error occurred in Panel B for a half-crack length 
of 2 inches. 

The calculated crack lengths and crack propagation rates shown 
on figures 16 and 17 are a function of debond growth. On figure 18 
the debond aspect ratio, b/a, was plotted against the half-crack 
length for the two panels. The symbols on the figure indicate values 
of the debond aspect ratio obtained experimentally (see figure 10). 

As evident on the figure, the calculations indicated that the debond 
aspect ratio increases rapidly, especially for Panel A, before the 
half-crack length reaches 0.2 inch. Hence, the debond grows before 
the crack does . 

Of the two panels. Panel A exhibits the largest debond growth, and 
at a half-crack length of 2 inches the predicted debond aspect ratio was 
determined experimentally. For Panel B, the analysis predicted a debond 
aspect ratio of 0.45 at a half-crack length of 2 inches. In contrast, 



crack propagation rate da/dN, inches/cycle 
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Fig. 17. Experimental and calculated crack propagation rates 



debond aspect ratio b/a 
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Fig. 18. Debond aspect ratio versus crack length 
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the debond aspect ratio obtained experimentally was about 0.14. The 
discrepancy between the calculated and experimental values may be 
linked to the magnitude of strain energy release rate as the debond 
propagates in Panel B. Values of G in Panel B ranged from 
0.63 in-lbs/in at a half-crack length of 1 inch to 0.53 in-lb/in at a 
half-crack length of 2 inches (in contrast G values for Panel A 
ranged from 6.0 to 2.45 in-lbs/in). The values of G for Panel B 
were below values of G used to determine the empirical constants in 
equation (5). A few exploratory tests of the type performed in Appendix 
A revealed that a threshold value of G may exist. Below this threshold 
value debonding does not occur. Although proving the existence of a 
threshold was beyond the scope of this dissertation, if it does 
exist the analysis of Panel B would predict a debond aspect ratio much 
closer to the experimental value. 

The comparison between calculated and experimental values of 
predicted crack length, crack propagation rate and debond aspect 
ratio shown on figures 16 through 18 showed that the analysis has 
potential to predict crack growth in metals reinforced with composite 
materials. However, a true assessment of the accuracy of the analysis 
can be made only after a more extensive data base is developed. 



CHAPTER VII 


PARAMETRIC STUDIES ON CRACK AND DEBOND GROWTH 

To give insight about crack and debond growth in reinforced 
systems, the analysis was performed for reinforced systems with 
several different adherend thicknesses, debond sizes, and crack 
lengths. The metal adherend thicknesses studied were 0.05, 0.10, and 
0.15 inch; the composite adherend thicknesses were 0.025, 0.05, and 
0.075; the half-crack lengths were 0.5, 1.0, 1,5 inches; and the aspect 
ratios of the debond areas were 0.001, 0.5, and 1.0. For the various 
combinations of adherend thicknesses, crack lengths, and debond 
aspect ratios, the stress intensities, strain energy release rates, and 
remote stress that caused nonlinear behavior of the adhesive were 
calculated and are shown in Table 2. 

The first two columns of the table give the metal and adherend 

thicknesses, t„ and t . The third and fourth columns give the 
me 

half-crack length, a, and the debond aspect ratio, b/a. The fifth 
column gives the stress intensity in the metal adherend normalized by 
the remote stress, k/s. The sixth column gives the strain energy 
release rate at the debond front (along the longitudinal centerline) 
normalized by the remote stress squared, 6/s^. The last column gives 
the remote stress that would cause the adhesive layer in the reinforced 
system to behave nonli nearly. 
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TABLE 2 

CALCULATED VALUES FOR PARAMETRIC STUDY 


t 

m 

(in) 

(^n) 

a 

(in) 

b/a 

k/s 

G/s^ 

(in‘*/lbs) 

^yield 

psi 

0.050 

1 

0.025 

1 

0.50 

0.001 

0.171 

12.4E-10 

19,600 





0.50 

0.50 

0.240 

8.2E-10 

25,200 





0.50 

1.00 

0.289 

4.4E-10 

25,600 





1.00 

0.001 

0.170 

12.6E-10 

17,600 





1.00 

0.50 

0.286 

9.5E-10 

23,900 





1.00 

1.00 

0.366 

4.9E-10 

23,600 





1.50 

0.001 

0.171 

12.4E-10 

16,500 





1.50 

0.50 

0.322 

9.5E-10 

23,800 



0.025 

1.50 

1.00 

0.426 

5. IE-10 

23,700 



0.050 

1 

0.50 

0.001 

0.138 

7.3E-10 

25,000 





0.50 

0.50 

0.190 

5.6E-10 

21 ,700 





0.50 

1.00 

0.230 

3.0E-10 

29,500 





1.00 

0.001 

0.136 

7.2E-10 

22,800 





1.00 

0.50 

0.221 

6.0E-10 

28,100 





1.00 

1.00 

0.286 

3.3E-10 

27,000 





1.50 

0.001 

0.137 

7.2E-10 

21 ,100 





1.50 

0.50 

0.264 

6.0E-10 

27,600 



1 

0.050 

1.50 

1.00 

0.329 

3.3E-10 

27,200 



0.075 

0.50 

0.001 

0.125 

5.7E-10 

28,200 

1 

0.050 

0.075 

0.50 

0.50 

0.170 

4.8E-10 

33,800 
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TABLE 2-Continued 


(in) 

(in) 

a 

(in) 

b/a 

k/s 

(in'') 

G 's^ 

(in**/lbs) 

^yield 

psi 

0.1 

050 

0.1 

075 

0.50 

1.00 

0.206 

2.4E-10 

31,900 





1.00 

0.001 

0.123 

5.6E-10 

26,500 





1.00 

0.50 

0.196 

4.7E-10 

30,600 



i 

1 


1.00 

1.00 

0.254 

2.6E-10 

29,100 





1.50 

0.001 

0.123 

5 6E-10 

23,800 





1.50 

0.50 

0.217 

4.8E-10 

30,000 

0.( 

150 

0.( 

175 

1.50 

1 .00 

0.291 

2.7E-10 

29,300 

0.100 

0.( 

125 

0.50 

0.001 

0.254 

41.5E-10 

11,200 





0.50 

0.50 

0.334 

24.6E-10 

15,700 





0.50 

1.00 

0.385 

11.4E-10 

17,600 





1.00 

0.001 

0.256 

46.1E-10 

9,800 





1.00 

0.50 

0.407 

28.9E-10 

14,200 





1 .00 

1.00 

0.498 

13.5E-10 

15,400 





1 .50 

0.001 

0.257 

46.2E-10 

9,100 





1.50 

0.50 

0.462 

29.9E-10 

13,700 



o.c 

)25 

1.50 

1.00 

0.583 

14.2E-10 

15,100 



0.050 

0.50 

0.001 

0.202 

23.9E-10 

14,600 





0.50 

0.50 

0.263 

16.2E-10 

18,900 





0.50 

1.00 

0.307 

8.1E-10 

19,900 


1 

1 



1.00 

0.001 

0.202 

25.1E-10 

13,100 

0.1 

00 

O.C 

)50 

1 .00 

0.50 

0.309 

18.5E-10 

17,400 
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TABLE 2-Continued 



0.100 


0.100 

0.150 


a b/a k/s G/s^ 

(in) (in) (in^) (mVlbs) psi 


0.050 1.00 1.00 

1.50 0.001 

1.50 0.50 

0.050 1.50 1.00 

0.075 0.50 0.001 

0.50 0.50 

0.50 1.00 

1.00 0.001 

1.00 0.50 

1.00 1.00 

1.50 0.001 

1.50 0.50 

0.075 1.50 1.00 

0.025 0.50 0.001 

0.50 0.50 

0.50 1.00 

1.00 0.001 

1.00 0.50 

1.00 1.00 

1.50 0.001 

1.50 0.50 

0.025 1.50 1.00 


0.385 

9.4E-10 

12,379 

0.201 

25.0E-10 

16,800 

0.345 

19.0E-10 

16,800 

0.444 

9.9E-10 

16,700 

0.178 

17.6E-10 

16,900 

0.230 

12.5E-10 

21 ,300 

0.269 

6.4E-10 

21 ,800 

0.177 

18.0E-10 

15,200 

0.266 

14.0E-10 

19,400 

0.333 

7.4E-10 

18,500 

0.176 

17.9E-10 

14,600 

0.294 

14.4E-10 

18,600 

0.382 

7.4E-10 

18,100 

0.313 

78.3E-10 

8,300 

0.397 

40.9E-10 

12,400 

0.446 

18.0E-10 

14,800 

0.323 

95.6E-10 

7,000 

0.493 

51 .3E-10 

10,900 

0.584 

22.3E-10 

12,500 

0.325 

99.0E-10 

6,500 

0. 565 

54.1E-10 

10,400 

0.689 

23.7E-10 

12,100 


0.150 
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TABLE 2-Continued 


(in) 

'c 

' in) 

a 

(in) 

b/a 

k/s 

G/$2 

(inVlbs) 

^yield 

psi 

0.' 

150 

0. 

050 

0.50 

0.001 

0.253 

46.8E-10 

10,600 





0.50 

0.50 

0.318 

28.9E-10 

14,400 





0.50 

1.00 

0.362 

13.9E-10 

16,200 





1.00 

0.001 

0.256 

52.9E-10 

9,400 





1.00 

0.50 

0.379 

35.2E-10 

12,900 

1 




1.00 

1.00 

0.460 

17.1E-10 

13,400 





1.50 

0.001 

0.256 

53.7E-10 

8,700 



i 

! 


1.50 

0.50 

0.425 

37.1E-10 

12,400 

0.1 

150 

0.( 

150 

1.50 

1.00 

0.534 

18.2E-10 

12,900 

0.150 

0.: 

75 

0.50 

0.001 

0.222 

34.4E-10 

12,300 





0.50 

0.50 

0.278 

22.5E-10 

16,100 





0.50 

1.00 

0.318 

11.3E-10 

17,600 





1.00 

0.001 

0.223 

37.3E-10 

11 ,100 





1.00 

0.50 

0.325 

26.9E-10 

14,500 





1.00 

1.00 

0.397 

13.7E-10 

14,400 





1.50 

0.001 

0.222 

37.6E-10 

10,300 





1.50 

0.50 

0.360 

28.2E-10 

14,000 

0.1 

50 

0.7 

'5 

1.50 

1.00 

0.457 

14.5E-10 

13,700 
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Stress Intensity Factor, Crack Growth Rate 

Figures 19, 20, and 21 show the effects of composite adherend 
thickness, debond aspect ratio, and crack lengths on stress intensity 
for the different thickness metal adherends. On the figures circles, 
diamonds and squares represent stress intensity ^/alues for 0.025, 0.05, 
and 0.075 inch thick composite adherends respectively. Open, half- 
shaded, and shaded symbols represent stress intensity values for 
debond aspect ratios of 0.0, 0.50, and 1,0 respectively. On the figures 
the stress intensity normalized by the remote stress was plotted against 
the half-crack length. 

For all three metal adherend thicknesses, the figures show 
consistent trends. If the debond aspect ratio is small, the stress 
intensity is not significantly affected by either the crack length or 
the thickness of the composite reinforcement. However, as the 
debond size increases, the stress intensities increase significantly 
with longer crack lengths and thinner composite reinforcement. 

The effects of debond size, composite adherend thickness, and 
crack length become more pronounced for the thicler metal adherends. 

Strain Energy Release Rate, Debond Propagation 

Figures 22, 23, and 24 show the effect of composite adherend 
thickness, debond aspect ratio, and crack length on the strain energy 
release rate at the debond front (along the longitudinal centerline 
of the reinforced system). As on figures 19, 20, and 21, different 




50 1.0 15 

a, inches 


ig. 19. Stress intensities for 0.050 in thick 
metal adherend 
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0.5 1.0 1.5 

a, inches 


Fig. 22. Strain energy release for 0 05 in thick 
metal adherend 
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Fig. 23. Strain energy release for 0.10 in thick metal 
adherend 
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Fig. 24. Strain energy release for 0.15 in thick 
metal adherend 
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symbols indicate different composite adherend thicknesses, and 
different amounts of shading indicate different aspect ratios. On 
the figures the strain energy release rate normalized by the remote 
stress squared was plotted against the half crack length. 

Examination of the figures revealed that for all metal adherend 
thicknesses the thickness of the composite had the most pronounced 
effect on the strain energy release rate. The thinner the composite 
adherend the higher the strain energy release rate and the more 
likely debonding would occur. The debond aspect ratio had the next 
most significant effect. The larger the debond aspect ratio the lower 
the strain energy release rate. Of all the parameters the crack 
length had the smallest effect on the strain energy release rate. 

As in the case for stress intensities, the effects of debond 
size, composite adherend thickness, and crack length are more 
pronounced for the thicker metal adherends. In fact on figure 24 
the energy release rates for 0.025 inch composite reinforcement with no 
debond reinforcing a 0.15 inch metal were so great for all crack 
lengths that the energy release rates were off scale in the figure. 
Evidently, the most severe case for debond occurs with a thick 
metal adherend reinforced with a thin composite sheet with no 
debonding between adherends. 

Nonlinear Effects 

Figures 19 through 24 were generated by assuming that the adhesive 
behaved linearly. Hence, because they were normalized with respect 
to the remote stress or its square, they can be used to estimate the 
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stress intensity and strain energy release rate for any remote stress 
that does not produce nonlinear behavior of the adhesive in the 
reinforced system. For the different parameters studied. Table 2 
gives the values for the remote stresses that cause nonlinear behavior 
of the adhesive. As evident from the table nonlinear behavior can 
occur at relatively low remote stresses. For value'; in the table, 
the lowest value of remote stress to cause nonlinear behavior occurred 
for the thickest (0.15 in) metal adherend with a 1.5 inch crack that 
was reinforced with the thinnest (0.025 in) composite adherend and no 
debond. For this reinforced system a remote stress of 6,000 psi 
caused nonlinear behavior of the adhesive. However, for practical 
purposes, the remote stress applied to this system could be as high 
as 52,000 psi. To investigate the effects of the nonlinear adhesive 
on the crack and debond growth predictions, the analysis was conducted 
using the preceding parameters for both a linear adhesive and nonlinear 
adhesive. 

For the linear analysis an effective shear modulus (see Appendix 
B) of 65,000 psi was used, while for the nonlinear analysis an effective 
shear modulus of 65,000 psi was used until the remote stress reached 
6,000 psi after which an effective shear modulus of 36,000 psi was 
used for nonlinear elements of region B (see page 38). With the use 
of the two different analyses, the stress intensity and the strain 
energy release were found as 
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adhesive 


G 

1 inear 

18,800 

26.1 

nonl inear 

16,700 

26.1 

percent error 

13 

0 


For this example, the nonlinear analysis predicted the stress 
intensity 13 percent below that predicted by the linear analysis. 
However, the predicted strain energy release rate was the same for 
both analyses. 

At first glance the linear analysis would then accurately 
predict the debond growth and conservatively estimate the crack growth. 
But by the use of equations (4) and (5) and the above values the debond 
and crack growth rates were found at 4.2 inches/cycle and 1.56E-04 
inches/cycle respectively. Hence, the debond propagates much faster 
than the crack. As the debond grows, the magnitude of the stress 
intensity from the linear and nonlinear analysis converges. In fact, 
because the debond grows so much faster than the crack, before the 
crack extends any appreciable amount the stress intensities from the 
two analyses predict the same crack growth rates. In addition, because 
the example exhibits the most significant nonlinearity of all the cases 
considered in Table 2, the linear analysis, and hence figures 19 
through 24 can be used to estimate crack and debond growth even when 
the adhesive behaves nonlinearly. 
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Prediction of Crack and Debond Growth 

To predict the crack and debond growth in a reinforced system, 
the computer code discussed in Appendix F should in general be employed. 
However, figures 19 through 24 can also be used to estimate both the 
debond and crack growth for a variety of reinforced systems in the 
following manner. First, for the given adherend thickness, crack 
length, debond size, and remote applied stress, the strain energy 
release rate, G, and the stress intensity, k^, can be estimated 

from figures 19 through 24. Then, equations (4) and (5), the rrack 
and debond growth equations, can be used to estimate the number of 
applied load cycles to extend the crack by, for example, 10 percent 
and to extend the debond length (along the longitudinal axis of the 
reinforced system) by 10 percent. The smallest value of the applied 
load cycles required to produce these extensions is used with equations 
(4) and (5) to predict the extension of crack and debond growth. Next, 
these extensions are added to the original crack and debond length 
and the entire process is repeated until the stress intensity reaches 
56,000 psi-in^ (fracture occurs) or the desired crack length or 
number of applied load cycles is reached. 



CONCLUSIONS 


The failure mode of cracked metal sheets that are reinforced 
with composite is crack propagation in the metal sheet. Analysis of 
the crack growth is complicated by the development of a debond near 
the crack. Herein, an analysis was developed to predict both the debond 
and crack growth in a reinforced system. The analysis was predicated 
on the use of strain energy release rate to correlate debond growth. 
Empirical constants required for the correlation were developed 
from simple bonded specimens. The correlating equation for the 
debond growth was then used in a stress analysis that was based on 
complex variable Green's functions which were developed herein for 
cracked, isotropic sheets and uncracked, orthotropic sheets. The 
stress analysis was used to calculate the inplane and interlaminar 
stresses, the stress intensity at the crack tip, and the strain 
energy release rate at the debond front. By the use of the analysis, 
an iterative solution was developed that used the stress intensity and 
the strain energy release rate to predict the crack and debond growth 
on a cycle- by- cycle basis. 

To verify the analysis, tests were conducted on two different 
reinforced panels which exhibited different amounts of debonding. 

For both panels the predicted crack growth was within the accuracy 
of crack growth prediction in unreinforced metal sheets. Hence, 
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the analysis appears accurate within the bounds of existing fracture 
mechanics concepts. 

The analysis was used in a parametric study of the effects of 
boron/epoxy composite reinforcement on crack propagation in aluminum 
sheets. The study showed that the aspect ratio of the debond area 
has a significant effect on the crack propagation in the aluminum 
sheet. For small debonds the crack propagation rate is reduced 
significantly, but these small debonds have a strong tendency to 
enlarge. Debond growth is most likely to occur m reinforced systems 
that have a cracked metal sheet that is reinforced with a relatively 
thin composite sheet. 

The analysis can be used to predict crack growth in reinforced 
systems. Hence, the analysis can be applied in developing methods 
to repair damaged metal structure and to increase lives and payloads 
of metal structures by selective reinforcement. 



APPENDIX A 


DETERMINATION OF DEBOND CONSTANTS 

As discussed in Chapter II, debonding can be predicted in a 
bonded system with an equation of the type 

nj 

db/dN = C 2 (G) (5) 

where G is the strain energy release rate and and are 

empirical constants. The objective of this appendix is to determine 
the empirical constants for the reinforced system used in the 
experimental portion of this dissertation. 

Specimen Fabrication 

To determine and n^, several test specimens with the 
configuration shown in figure A.l were fabricated. The specimens 
consisted of 1-inch wide strips of 0.188 inch thick 2024-T3 aluminum 
bonded to 0.03 inch thick unidirectional boron/epoxy The strips were 
bonded with Shell EA-934 room curing adhesive. To maintain a constant 
adhesive thickness in the bond, 2 percent by volume of 0.004 inch 
diameter glass beads were added to the adhesive prior to bonding. 

The process used to bond the aluminum to the composite was as 
follows 
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Fig. A.l. Debond specimen configuration 
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SURFACE PREPARATION 
Aluminum 2024-T3 

1. Vapor degrease - perchloroethylene condensing vapors for 

5 to 10 minutes. 

2. Grit blast with 220 grit at 90 psig. 

3. Alkaline degrease - Oaklite 164 solution (9 - 11 ounces/ 

gallon of water) at 190 ±10° F for 15 minutes. Rinse 
immediately in large quantities of cold running water. 

4. Acid etch - place panels in the following solution for 

10 minutes at 150° ±5° F. 

Distilled water 30 parts 

Sulfuric acid (cone) 10 parts 
Sodium Dichromate 1 part 

5. Rinse - rinse panels in clear, deionized running water. 

6. Dry - air dry 15 minutes; force dry 10 minutes at 150° F 

±10° F. 

Boron/epoxy 

1. Vapor degrease as above. 

2. Grit blast with 220 grit at 30 psig. 

BONDING 

1. Bond within 4 hours of surface preparation. 

2. Coat surfaces of both adherends prior to bonding. 

3. Cure at room temperature under 15 psig ±2 psig pressure. 

4. Record date of bonding. 
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The method of surface preparation for the aluminum was taken from 
Cagle (1973) while the surface preparation for the boron/epoxy and the 
bonding method was developed by the author. The bonding process was 
verified with lap-shear strength tests of the bonded system. 

As shown in figure A.l photoelastic material was bonded to the 
surface of the boron/epoxy. The photoelastic material enabled tracking 
of the debond front in the fatigue tests of the specimens. 

Fatigue Tests 

The fabricated specimens were tested in a servo-hydraulic fatigue 
machine with a maximum load capacity of 10,000 pounds. All of the 
tests were conducted at a loading frequency of 10 Hz with a ratio of 
minimum to maximum load in the load cycle of R = 0.05. Duplicated 
tests were conducted for maximum loads of 5,000, 4,000, and 3,000 
pounds. 

As the specimens were tested, a debond developed at the change 
of cross section and propagated between the aluminum and boron/epoxy 
adherends. Throughout the tests the location of the debond front 
was indicated by an isochromatic that was observed by viewing the 
photoelastic material through a polarizing material. The location 
of the debond front is plotted against the number of applied load 
cycles on figure A. 2 for all of the tests. The results of these 
tests will be used with a stress analysis to determine the empirical 
constants c^ and nj. 
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1 OOK ZOOK 

applied load cycles, N 

Fig. A. 2. Debonding behavior 



Stress Analysis 


Because of the change of cross section of the test specimen, 
as an axial load, P, is applied to it, it bends. As shown by 
Timoshenko (1961) the equilibrium equation for a beam that exhibits 
both axial and bending deformation can be written as 


dM dw 

V = — + P— (A.l) 

dy dy 

where V is the shear, M is the bending moment, and P is the 
tensile axial load. To use equation (A.l) both the moment M and 
the axial load P were related to the deflection in the z-direction, 
with two equations given by Calcote (1969) as 


P 

dv 

d^w 


A 

A — 

1 1 j 

dy 

dy"" 

(A. 2) 


dv 

d^w 


M = 

dy 

dy^ 

(A.3) 


where A is the cross sectional area of the beam and dv/dy is the 
axial strain of the beam and 
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1 - (v V ) . 
' xy yx'i 


where i indicates the layer of a beam with k layers and z is 
measured from any reference surface. 

Eliminating dv/dy from equations (A. 2) and (A. 3) yielded an 
expression for the moment as 



d^w 

D 

11 dy2 


(A.4) 


Substituting equation (A.4) into equation (A.l) and differentiating 
once with respect to y, yielded a governing equation as 


^z 



d“w d^w 

Dll ^ + P 

dy** dy^ 


(A. 5) 


where q^ is a transverse distributed load acting on the beam. 

The boundary conditions for equation (A. 5) were determined 
from the end conditions of the test specimen installed in the test 
machine as shown schematically on figure A. 3a. With the assumption 




(b) beam model of specimen 


Fig. A. 3. Beam model for stress analysis 
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that the specimen is effectively fixed at both ends by the relatively 
massive test machine grips the boundary conditions were determined as 


w(-d) = 0 w{L) = 0 


dw ( - d ) dw ( L ) 

= 0 = 0 


dy 


dy 


(A. 6) 


At the change in cross section of the test specimen a local 
moment given as 


M = Pe^ (A. 7) 

where e^ is the distance between centroids of the two cross sections, 
was produced. With this moment the beam can be modeled as a symmetric 
beam with two different cross sections acted upon by a local moment 
M at the change in cross section and an axial load P as shown in 
figure A. 3b. For the model shown in figure A. 3b, equation (A. 5) 
was solved with conventional finite difference techniques (Ames 1971). 

To verify the analysis, the strains were determined by both the 
finite difference analysis and by experiments. The test specimen 
shown in figure A. 3a was analyzed. In cross section 2 the specimen 
was composed of four layers: the metal core, the adhesive layer, the 

composite cover, and the photoelastic material. The thickness, moduli, 
and Poison's ratio for each of these layers are given as 
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layer 

material 

thickness 

(in) 

elastic 

modulus 

psi 

Poisson’s 

ratio 

1 

aluminum 

0.188 

10 ^ 

0.30 

2 

adhesive 

0.004 

600,000 

0.40 

3 

composite 

0.030 

3 X 10 ^ 

0.21 

4 

photoelastic 

0.083 

390,000 

0.36 


With the preceding values and the surface of the metal with strain 
gages (see figure A. 4) as a reference plane Bjj, and Djj 

were found for section 1 (see figure A. 3) as 


A = 2.1 X 10^ Ibs/in B = 200,000 lbs 
11 11 

D = 2,500 Ibs/in 
11 

and for section 2 as 

A^j = 3 X 10® Ibs/in B^^ = 392,000 lbs 

Djj = 65,000 Ibs/in 

With the use of the previous values and an applied load P of 5,000 
pounds as shown on figure A. 4, the finite difference method was used 
to calculate deflections and curvatures of the beam. With the curva- 
tures the strains in the beam were calculated on the surface of the 
metal with (Cal cote 1969) 
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On figure A. 4 the solid line indicates calculated strains on the 
surface of the metal adherend while the circles indicate strains 
obtained experimentally with strain gages. As evident from the figure, 
the comparison is very good. Consequently, the developed stress 
analysis is adequate and can be used to calculate the strain energy 
release rates as the test specimens debond. 

Calculation of Strain Energy Release Rates 
The strain energy release rate can be calculated as 

lim r AW AU') 

G = S f (A. 8) 

Ab ^ 0 / Ab Ab j 

where W is the change in work done on the system, U is the change 
in internal strain energy, and Ab is a small extension of the 
debond. With the assumption that the debond extends at the maximum 
applied load, the work done as the debond extends can be calculated as 


AW 


^^•^axial ^ 


A6 


bending 


, + 

. 


MA00 


(A. 9) 


where A0G is the rotation of the beam at the debond front. The 
change in internal energy can be calculated as 


axial 


+ AU 


AU = AU 


bending 


(A. 10) 
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Substituting equations (A. 9) and (A. 10) into equation (A. 8) and 
rearranging slightly yields the strain energy release rate as 


[ '^axial 

G = 1 im < P 

Ab ^ 0 ) Ab 


AU 


axial 

Ab 


^bend 

+ P 

Ab 


MAOG 


bend 


Ab 


AU 


bend 


Ab 


(A. 11) 


The first two terms, which neglect bending, of equation (A. 11) can be 
calculated as shown by Roderick et al (1975) as^ (for a unit width) 


1 im 

Ab ^ 0 


r <5 

axial 

< p 

Ab 


AU , I t E P^ 

axial ( c c 

Ab \ 2E_r (t,,,E., + l_E.) 
j m m 111 til c c 


(A 12) 

The last three items, which are due to bending, can be calculated with 
the finite difference results in the following manner 

First, the axial deflection due to bending is calculated as 
shown by Den Hartog (1952) as 



The flexible adhesive and photoelastic material are neglected 
because they have little effect on the strain energy due to the axial 
load 


A 



no 


where the slope dw/dy was calculated at each nodal point in the 
finite difference approximation. Between the nodes the slope 
was assumed to vary linearly. The integration of equation (A. 13) 
was done piecewise over the length of the beam. The work done 
by the moment was expressed as 


M 



(A. 14) 


The strain energy due to bending was calculated as 


A 


1 1 


U = / 

^bend -d 


2(A D - B ) 
'■11 11 n ' 


dy 


(A. 15) 


where again the integration was performed in a piecewise fashion. 

Equations (A. 13) through (A. 15) were evaluated before and after 
a debond extension Ab. The quantities were subtracted to calculate 
the last terms of equation (A. 11). The result was added to equation 
(A. 12) to give the strain energy release at the debond front. In 
figure A 5 the solid lines show the calculated strain energy release 
rate plottedaqainst the debond lengths for values of applied loads of 
5,000, 4,000, and 3,000 pounds. The dashed line on the figure shows the 
strain energy release rate--negl ecting bendi ng--cal culated with 
equation (A 12) As evident from the figure, when the debond length 
IS greater than 0 5 inch or less than 4.5 inches, the contribution 
of bending to the strain energy release rate is small. Particularly, 



strain energy release rate, in-lb/in 
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Fig. A. 5. Strain energy release rate versus debond front 
location 
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when the debond length is about 2,5 inches long, the bending contribution 
2 

IS zero. Consequently, for a debond length of 2 5 inches, equation 
(A. 12) can be used to calculate the strain energy release rate 

Curve Fit for Empirical Constants 

The debond rates for the configuration analyzed were determined 
by taking the slope of the curves at a 2.5-inch debond length from 
figure A. 2. These experimental rates and the corresponding values 
of the calculated strain energy release rates are shown in Table A.l. 

The values in Table A.l were used to determine the empirical constants 
C 2 and n^ in equation (5) with a least squares fit (Wylie 1966). 

To perform the fit, equation (5) was written as 


log(db/dN) = log(c 2 ) + n^logCo) 


(A. 16) 


As a result of the curve fit, Cj and nz were found as 


c = 3.158 X 10 
2 


n - 3.616 

2 


Figure A. 6 shows the data as points and the fitted equation (5) 
as a solid line. As evident from the figure, the equation fits the 


2 

Actually, there are also two other deoond lengths where the 
strain energy release rate is zero, but they are located closer to 
the ends of the specimen where the analysis may be more inaccurate 
than at the 2.5-inch debond length. 



113 


TABLE A.l 

STRAIN ENERGY RELEASE RATES 
AND DEBOND PROPAGATION RATES 


load 

strain energy 

debond rate 


release rate, G 

db/dN 

lbs 

in- lbs/ in 

in/cycle 


5,000 

5.000 

4.000 

4.000 

3.000 
3,000 


2.15 

2.15 

1.38 

1.38 

0.77 


5.60 X 10 ^ 
5.20 X 10 “^ 
1.04 X 10 ”^ 
8.00 X 10 "^ 
1.26 X 10 "^ 
1.33 X 10 ’^ 


0.77 




debonding rate, in/cycl 
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Fig. A. 6. Correlation of debonding rates with strain energy 
release rates 
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data well. With the determined empirical constants, equation (5) 
can be used to predict debond propagation rates whenever the strain 
energy release rate can be determined. 



APPENDIX B 


ADHESIVE SHEAR DEFORMATION ASSUMPTION 

As mentioned in Chapter IV, the complexity of the analysis 
for a reinforced system is significantly reduced by assuming that the 
adhesive only undergoes shear deformation while the adherends only 
undergo extensional deformations that do not vary through the 
adherend thickness. Herein, the validity of these assumptions were 
investigated by comparing a one-dimensional solution that was based 
on these assumptions with a more rigorous two-dimensional finite 
element solution. Before the one-dimensional and two-dimensional 
solutions were compared the bulk properties of the adhesive were 
determined. 


Adhesive Bulk Properties 

To determine the bulk properties of the adhesive an appropriate 
test specimen was designed and fabricated in several steps. First, 
a female plastic mold was made from the specimen shown in figure B.la. 
Then, the adhesive liquid base and hardener were combined and cast 
into the mold. Next, after curing 24 hours in the mold the adhesive 
specimen was removed and cured an additional 5 days before it was 
handled. Finally, x-rays were taken of the specimen to locate air 
bubbles developed in the molding process. Specimens that contained 
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large air bubbles in critical areas were scrapped Six test 
specimens were fabricated in this manner but only three were acceptable 
for testing. 

Prior to the testing the specimens were instrumented with strain 
gages and linear variable differential transformers (LVDT) as shown 
in figure B.lb. The strain gages were used to obtain the longitudinal 
and transverse strain while the LVDT's were used to check the strain 
gages (the LVDT's were wired to eliminate bending effects in their 
readings while the strain gages would include bending in their 
readings). Discrepancies between the readings would indicate bending 
due to poor specimen alignment in the test machine Each instrumented 
specimen was placed in a servo-hydraulic test machine with a loading 
range of 2,000 pounds and a sensitivity of ±10 pounds. Then, each 
specimen was loaded to failure at a rate of 80 Ibs/sec. 

In this manner, the three specimens were tested and the results 
were nearly identical for all three of the tests. For all specimens, 
the longitudinal strain calculated from LVDT's agreed within 1.5 percent 
of the longitudinal strain gage reading and indicated the specimens 
were aligned properly in the test machine. The data obtained from the 
strain gages are shown on figure B.2 in the form of a stress-strain 
plot. On the figure the heavy solid line indicates the strain in the 
loading direction and the dashed line indicates the strain 
transverse to the loading direction e . As indicated by the light 

A 

solid lines on the figure, the stress-strain curve can be approximated 
by a bilinear stress-strain curve with a change of slope occurring 
at 4,200 psi. 



stress, ksi 



Fig. B.2. Adhesive stress-strain curve 
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In the initial linear region the following values were obtained 
from the curve: 

AOy = 3,000 psi, Ae^ = 0.005, = 0.002 

Using these parameters, the Poisson's ratio and the elastic modulus 
v;ere calculated as 


yx 


= 0.40 


Ao. 


Ae. 


= 600,000 psi 


With the assumption that the adhesive is isotropic the shear modulus 
was calculated as 


Q = = 215,000 psi 

2 (1 + V ) 

yx' 

Along similar lines, the material parameters in the second linear 
region were determined as 

V = 0.28 E = 190,000 psi G = 74,000 psi 

y 

The fracture of the specimen occurred at a stress of 6,600 
psi at an axial strain, c^, of 0.0215. With the preceding material 
property values for the adhesive, the one-dimensional and two- 
dimensional solutions were compared. 
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As a test case, the shearing stresses in the adhesive layer 
of the specimen shown on figure B.3a (the specimen is symmetric about 
the x-y plane) was calculated with both types of solutions for a plane 
stress state. Although the solution was for plane stress, for the 
self-equilibrating load system shown on figure A. la the stress 
distributions are identical for both plane stress and strain states 
(Timoshenko (1951)). Consequently, although the tests case considered 
a state of plane stress, the results are also applicable to a state of 
plane strain which may be more appropriate for a section taken through 
the thickness of the reinforced system shown on figure 9. 

One-Dimensional Solution 

Figure B.3b shows a freebody of the specimen shown in figure B.3a. 
On the figure P is half the load in a composite adherend, F(y) 
and q(y) are the load in the metal adherend and the shear flow in 
the adhesive at any point y, and t^^, 2t^, and t^^ are the 

thicknesses of the metal, composite, and adhesive respectively. 

The change in the load F(y) with respect to y is the shear 
flow in the adhesive layer given as 

dF(y) 

q(y) = (B.i) 

dy 

The shear stress in the adhesive is this shear flow divided by the 
width, w, of the specimen, given as 

q 1 dF(y) 

T(y) = — = (B.2) 

w w dy 



adhesive 


L composite 

(b) Freebody for one-dimensional solution 


Fig. B.3. Specimen configuration and freebody for 
one-dimensional solution 
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With the assumption that the adhesive only undergoes shear deformation 
and that this deformation does not change through the adhesive thickness, 
the shearing stress in the adhesive can be related to the shearing 
strain in the adhesive by the constitutive equation as 

T(y) = Y(y)G 3 j (B.3) 

The shearing strain in the adhesive can be related to the extensional 
displacements in the metal, ^nd the composite, u^(y), 

as 

V Jy) - v^(y) 

Y(y) = (B.4) 

^ad 

Substituting equation (B.4) into equation (B.3) and that result into 
equation (B. 2) yields 

dF(y) 

= \y^{y) - V (y)] (B.5) 

Equation (B.5) can be differentiated with respect to y to yield 


d^F 

dy^ 


dv.Jy) dv.fy) 



(B.6) 
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when the derivatives of the extensional displacements and 
with respect to y are the extensional strains in the metal and 
composite adherends respectively. These strains were related to the 
extensional loads in the adherends by 

dv^(y) F{y) 

e = = 

^ dy wE t 

^ mm 


dv^(y) P - F(y) 

£ = = (B.7) 

dy wE^t_ 

c c 


where E^^ and E^ are the moduli of the metal and composite respec- 
tively. Substituting equations (B.7) into equation (B.6) yielded a 
second order, nonhomogeneous differential equation in F(y) as 

d^F(y) 

- aF(y) = B (B.8) 

dy 


where 


and 


"ad 


r 


a = 


ad i 


^m^in 


1 1 


^c^cj 


-PG 


ad 


-sG 


ad 


B = 


wt, .E^t^ 
ad c c 


'ad^c 
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where s is the stress in the composite and is the shear modulus 

of the composite. 

Solving equation (B,8) yielded a complete solution as 


t/ay -vSy 3 

F(y) = Cje + CjO + - 

a 


(B.9) 


With referral to figure B.3a, the boundary conditions for equation (B.9) 
are 


y = 0 F(y) = F(0) = 0 


y = oo F(y) IS bounded 


(B.IO) 


With the use of these boundary conditions, the constants in equation 
(B.9) were found to be 


S = 0 


3 

"2 = a 


Thus, equation (B.9) becomes 


swG 


F(y) = 


ad 


at .E 
ad c 



(B.ll) 


By the use of equation (B.2), the shear stress in the adhesive was 
calculated from equation (B.ll) as 



Finite Element Solution 


A finite element computer program, PLANE, was used to calculate 
the shearing stresses in the adhesive layer in the specimen shown in 
figure B3.a. PLANE is an elastic-plastic, two-dimensional program 
which uses constant strain and linear strain triangular elements. 

PLANE was developed by the Grumman Aircraft Corporation for the 
National Aeronautics and Space Administration and was documented by 
Armen and Levy (1962). The mesh used for the analysis is shown on 
figure B.4 and contains 1,522 degrees of freedom. The triangles 
are predominately linear strain triangles which allow linear variations 
in the stresses and strains through the elements. Each of the adherends 
is modeled with several elements through the thickness thus allowing 
variations of extensional and shearing stresses through the thickness. 

In contrast, the one-dimensional analysis assumed uniform extensional 
stresses and no shearing stresses through the thickness of each 
adherend. The adhesive layer was modeled by one layer of elements 
which allowed linear variation in stresses through the adhesive 


thickness. 
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One-Dimensional Versus Finite Element Solution 

To compare the results of the two solutions, the specimen 
configuration shown in figure B.3a with a width of w = 1.0 inch 
and the following parameters 



metal 


adhesive 

composite 

thickness 

t = 0.1 in 
m 

^ad 

= 0.004 in 

t^ = 0.1 in 

modul i 

E_ = 10.3 X 10® psi 
m 

^ad 

= 215,000 psi 

E^ = 30 X lo' 


and an applied stress s of 1,000 psi was analyzed with both solutions. 
The parameters used in these solution were typical for constituents 
used in the reinforced system discussed in Chapter II. The adhesive 
shear stress calculated from the one-dimensional solution equation 
(B.12) was plotted as a solid line against y (distance from the 
edge of the metal adherend) on figure B.5. On the same figure the 
circles indicate values of the adhesive shear stress calculated 
from the finite element solution. As evident from the figure the 
one-dimensional solution gives shear stress values twice as high as 
those obtained from the finite element solution near the edge of the 
metal adherend (y = 0). Evidently, the shearing deformation of the 
adherends which was not accounted for in the one-dimensional solution 
has a significant effect on the values of the shear stresses. 

To account for the adherend shear deformation and still use the 
simplified one-dimensional analysis, an effective shear modulus 



shear-stress in adhesive 
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Fig. B.5. Finite element versus one-dimensional solution 
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was introduced. The magnitude of was determined by equating 

the maximum shear stresses in the adhesive calculated by the finite 
element solution to the expression for the maximum shear stress from 
the one-dimensional solution from equation (B.12) as 


T 

max 

Finite Element 




(B.13) 


Solving equation (B.13) for G yielded an expression for an 
effective shear modulus Gg^^ as 



With this value for G in equation (B.12) yields a corrected one- 
dimensional solution for the adhesive shear stress as 


T(y) 


sGeffWe- 


^ad^c 


(B.15) 


For the sample analysis, equation (B.14) yielded G^^^ as 64,000 
psi. Using this value for equation (B.15) was plotted against 

y as a dashed line on figure B.5. The agreement between equation (B.15) 
and finite element results indicated by the circles on the figure is 


excellent. Consequently, the assumptions made in Chapter IV, that the 
adherends only undergo extensional deformation while the adhesive 
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only undergoes shear deformation, can be used to accurately predict 
shear stresses in the adhesive with a simplified analysis if an 
effective shear modulus for the adhesive is used in the calculations. 

For analysis of the reinforced system shown in figure 2 the 
effective modulus was determined for a range of adherend thicknesses 
for the adhesive used in the reinforced system both for before 
and after the adhesive yields. To determine the values of effective 
shear moduli numerous finite element solutions were run with different 
adherend thicknesses for both the initial bulk adhesive shear modulus 
of 215,000 psi and the bulk shear modulus after yielding of 74,000 
psi. The results of these calculations are given in table B.l. In 
the table the maximum shear stress calculated from the finite element 
solution and the effective shear modulus are tabulated for the 
different adhesive thicknesses. 

As shown in table B.l the value of the effective modulus for the 
initial shear modulus of the adhesive does not vary much with adherend 
thicknesses. The average value for is about 65,000 psi and 

is within 3 percent of any of the calculated values. 

Also, as shown in table B.l, the value of the effective modulus 
for shear modulus of the adhesive after yielding also has little 
variation with adherend thicknesses. The average value in this case 
is about 36,000 psi and is within 3.3 percent of any of the calculated 
values. 

As evident from the previous discussion the reinforced system 
can be analyzed by assuming that the adherends only undergo extensional 



132 


deformation while the adhesive only undergoes shearing deformation 
if an effective shear modulus of 65,000 psi is used for the adhesive 
before the adhesive yields and an effective shear modulus of 36,000 psi 
is used for the adhesive after yielding. 



TABLE B.l 

EFFECTIVE SHEAR MODULI FOR REINFORCED SYSTEM 
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APPENDIX C 


REMOTE STRESSES IN THE ADHERENDS 

For equilibrium to exist, the macroscopic stresses applied to the 
reinforced system must be balanced by the stresses in the adherends 
of the system. This relationship can be written in vector form as 


• 1 


am 

XX 


ac 

XX 



A = 

om 

yy 

A + 
m 

""yz 


^xy 

system 

am 

xy 

metal 

ac 

xy 

compos i te 


(C.l) 


where A represents the area of the reinforced system and A^^ and 
Aj. represent the area of each element. For the case where stress 
was applied in only the y-direction. 


S =0 


Sy = applied load/A 


s =0 

xy 


(C.2) 


With strains uniform through the thickness, (C.l) can be rewritten for 
a unit width as 
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where the C and D are the stiffness matrices for an isotropic 
and orthotropic material in plane stress, respectively, as given by 
Zienkiewicz (1971) as 





136 


with E, G , and \) denoting the extensional modulus, the 
yx yx 

shear modulus, and Poisson's coefficient respectively. With simple 
matrix algebra, equation (C.3) can be solved and the strains determined. 
These strains in turn can be used with the appropriated stiffness 
matrices to calculate the stresses in the metal and composite 
adherends as 




(C.7) 



APPENDIX D 


GREEN'S FUNCTION FOR THE CRACKED SHEET 

The solutions for the interlaminar stresses and the strain energy 
release at the debond front developed in Chapters IV and V respectively 
require Green's functions for both displacements and stresses. As 
pointed out by Dennemeyer (1968), the solution for a concentrated 
load acting on a body can be used as a Green's function. 

Herein, the solution for four concentrated loads that are 
symmetric with respect to a crack in an isotropic sheet was developed. 

The solution was based on elasticity theory using complex variable 

theory as developed by Muskhelishvili (1975). The solution was 
predicted on the assumptions that the cracked sheet is infinite 
isotropic, and can be described by a plane stress or strain analysis. 

As shown by Muskhelishvili (1975), both the stresses and displacements 
in a cracked sheet can be expressed in terms of two stress functions, 
^(z) and z), as 

= Real{$(z) + $(z) - [z$'(z) + 'i'(z)]} (D.l) 

Oy = Real{$(z) + $(z) + z’$'(z) + 'i'(z)} (D.2) 
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a = Imag{z$' (z) + 4 '^z)} 
xy 


(D.3) 


2G(u + iv) = n<|)(z) - z$‘ (z) - ii;(z) 


(D.4) 


where a , a and a are stresses, u and v are displacements, 
X y xy 

the bar over the functions denotesthe complex conjugate, i is the 
square root of -1 and 


d<j)(z) 

$(z) = <})' (z) = 

dz 


di]j{z) 

ftz) = ip' (z) = 

dz 


n = 3 - 4v plane strain 
3 - V 

n = plane stress 

1 + V 

is the Poisson's ratio. 

The stress functions, $(z) and f(z), for the cracked sheet 
under four concentrated loads as shown on figure D.lc were constructed 
by superimposing the stress functions for a cracked sheet under two 
different loading conditions. The first condition which is shown 
on figure D.la has concentrated loads acting on the cracked sheet plus 
a stress distribution applied to the crack surface. This stress 
distribution was equal to the stress distribution which exists along 
the x-axis for an uncracked sheet with concentrated loads. This stress 






j to formulate 6 


tions for a cracked sheet 


reen's func 
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distribution closes the crack and, therefore, effectively eliminates 
it. As a result, figure D.la represents an uncracked sheet with 
concentrated loads acting on it. The second condition which is shown 
on figure D.lb is for a stress distribution acting on the surface of a 
crack. This stress distribution has the same magnitude, but is 
of opposite sign to the stress distribtuion applied to the crack in 
the first condition. The development of the stress functions for 
these two conditions follows. 

The stress functions for figure D.la were developed by super- 
imposing the solution for single concentrated forces that act in 
different quadrants. Fora point, z^, in the first quadrant the 
solution was given by Muskhelishvili (1975) as 


1 


$(z) = 



H'(z) = Sn - 
z 



(D.5) 


where 


X + lY 

S = 

2it( 1 + n)t|^ 

and X and Y are the load components in the x and y directions 
respectively, it = 3.1459, and t^^ is the thickness of the sheet. 

Equations (D.5) were generated for the second, third and fourth 
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quadrants by replacing S and with -S and -z^, -S and -z^, 
and S and z^ respectively. These stress functions were then 
superimposed to form the functions 


$(z) = S 


'Vfz) = S<n 


+ S<n< 


1 

1 • 

. 


" ■ "o. 

( ’ 

’ ) 



• 1 

’ 1 




> + S 


1 1 


'z + z - z„ , 
0 0 , 


(D.6) 


-Z, 


(Z - Z^)=^ (z + z,)^ 


-z. 


' . 1 ! 




(D.7) 


where upon differentiation equation (D 6) becomes 


1>'(z) = -S 


(z * 


(Z - Z„)^ 


ys 


(z-z^)^ 


(D.8) 


and upon integration equations (D.6) and (D.7) become 
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<|)(z) = S [Log(z + Zq) - Log(z - z^) ] 


+ S [Log(z + z^) - Log(z - z^) ] 


(D.9) 



(D.IO) 

Equations (D.6) through (D.IO) are the required stress functions, 
derivative, and integrals to compute stresses and displacements from 
equations (D.l) through (D.14) for figure D.la. 

The stress functions for load condition 2 shown on figure D.lb 
were obtained in the following manner. Following Muskhelishvili (1975), 
a new stress function Q{z) was introduced which is related to the 
previously discussed stress function by the equation 


Tlz) = q(z) - $(z) - z4>* (z) 


(D.ll) 


where 


n(z) = n(z) 
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The required stress functions, $(z) and fi(z), can be 
determined for pressure acting on a single crack as (Muskhelishvili 
1975) 


1 


%{z) /" 

^ 2T7iI(z) -a 


I(t)p{t)dt 


t - z 


q(t)dt 


2iTi -a t - z 


(D.12) 


1 


fl (z) = 

^ 27:il(z) -a 


I(t)P(t)dt 


t - z 


q(t)dt 


2iTi -a t - z 


(D.13) 


with 


crack length 

I(z) = /z^ - a^ a = 

2 


and 


P(t) 


1 

2 




q(t) 





(D.14) 


Oy and T^y are normal and shearing stresses acting on the crack 
surfaces respectively (the plus sign indicates the upper surface of 
the crack while the minus sign indicates the lower surface.) P(t) 
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and q(t) are total normal and shearing pressures which act on the 
crack surfaces. As mentioned previously, the stresses which acts 
on the crack surface in figure D.lb is equal to the magnitude but of 
opposite sign to the stress caused by four point loads acting on an 
uncracked sheet. Therefore, the functions, $(z) and 4'(z), shown 
in equations (D.6) and (D.7), which are solutions for point loads on a 
continuous sheet, can be used to find P(t) and q(t). 

The normal and shearing stresses acting on the crack were 
calculated (Muskhelishvili 1975) in terms of equations (D.6), (D.7), and 
(D.8), as 


cJy - iTj^y = <I>(z) + $(z) + z$’(z)" + V(z) (D.15) 

Substituting equations (D.6), (D.7), and (D.8) into equation (D.15) 
yields 




oi(z,Zq) + a{z,z^) 


(D.16) 


or 


Oy - iT^y = 2Real (a(z,z^)) 


(D.17) 


where 


-4z, 


a(z,Zg) = S 


z2 - z. 






0 ) 
(D.18) 


2 


+ 


+ 
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Equation (D.16) shows that no imaginary term exists. Consequently, 

T is zero. Also, a has the same value along the x-axis independent 
xy j 

of the direction in which the x-axis is approached. Consequently, 

= Oy”. Therefore, with consideration given to the above statements 
and equations (D.14), the total pressures on the crack are 


P(t) = + aTtTi^ 

q(t) = 0 


(D.19) 


As a result of the above simplifications, equations (D.12) and (D.13) 
become equal to each other and are expressed as 


1 

4>(z) = Q{z) = 

27ril(z) 


I(t)P(t) 

/ a dt 

-a t - z 


(D.20) 


The integral in equation (D.20) was evaluated by contour 
integration along the contour shown in figure D.2 by using the 
Residue theorem given as (Wylie 1960) 


f f(z)dz = 2tTi E Residues 
c 

where the residue for simple poles is given as 



(D.21) 


Residue = lim 

(M - 1) ! z -*■ a 


M 

(z - a) f(z) 


(D.22) 
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Fig. D.2. Path for contour integration 
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The contour shown in figure D.2 can be broken into several sections. 
Thus, equation (D.19) becomes 

/ A(0 + / A(0 + / A(c) - / A(c) + / A(c) 

r, ^ 

- / A(?) + f A(c) = 2TTi E Res (D.23) 


where 

-a^ P(5) 

A(c) = d^, c = t + iS 

K - z 


On figure D.2 the integral is evaluated along the contour as R 
approaches infinity and e approaches zero. Each of the contours 
can be expressed as follows. 


/ A(?) = lim 

r, . . 


lU) P(c) 

/ dc 

C - z 


Let ? + a = ee 


i0 


dc = iee^®d0 


ie \/e"e^® 


2aee^® P(ce^® - a)e^®d0 


•• / A(c) = lim / 


e -»■ 0 


21T 


ee^® - a - z 


(D.24) 
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likewise for r, let ? - a = ee 


i0 


ie^e^e^® + 2aee^® P(ee^® + a)e^®d0 


.. , .-TT 

A(c) = Km \ 


«’® + a - z 


£ -*■ 0 


(D.25) 


for F 3 


I(t) P{t)dt 

f A(0 = f ® 

Fa -a t - z 


(D.26) 


likewise 


I(t) P(t)dt I(t) P(t) 

/ A(d = / ® = - / ® dt 

-a t - z -a t - z 


(D.27) 


- + 

Noting that Muskhelishvili (1975) showed that I{t) = -I(t) for 

^ < a, equations (D.26) and (D.27) can be combined to give the 

integral on the left hand side of equation (D.20) as 

+ 

I(t) P(t)dt 


/ A(^) = / A(c) = 2/a (D.28) 

Fj F^ -a t - z 



149 


Noting that 


I(t) 


= I(t) for c > a 


the sum of the next two segments of the contour integral equals 
zero, because 


/ A(c) = I 


R 


I(t) P(t)dt 


I(t) P(t)dt 


/ A(c) = 

t-z r R t-z 

6 


so that 


/ A(c) + / A(0 = 0 


(D.29) 


The last contour segment can be evaluated as 

K?) P(c)d? 

/ A(c) = lim / 

^7 R 00 C ■ z 


Let c = Re^®,d^ = Rie^®d0 


" / A(^) = lim / 

^7 R ^ 00 


21T 


Ri \fRV^®~~I^P(Re''®)d0 


Re^® - z 


(D.30) 
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Substituting equations (D.24), (D.25), (D.28), (D.29), and (D.30) 
into equation (D.23) and rearranging the terms yield the desired 
integral as 


I(t) P(t)dt 


1 


ie\/e^e^^“ 2aee^®P(ee^® - a)dG 


= Tfi E Residues - lim - / 

t - z n ? 

e ^ 0 ^ 


ee^® - a - z 


1 

lim - / 

n 2 TT 

e ^ 0 


-TT 


ie^e^e^Q + 2aee^® P(ee^® - a)e^®d 0 


ee^® + a - z 


1 

- lim - / 

2 ® 
R ^ 00 ^ 


21T 


Ri^^R^e^^® - a 2 P(Re^®)d0 


Re^® - z 


(D.31) 


where P(t) is given by equation (D.19). As an example, equation 
(D.31) was calculated for the first term of P(t). For this case 
the expression to consider was chosen as 


1(C) P(c) -41(c) Sz^ -4 SZq 

L_ = , p(t) = 

C - z (c== - - z) ^ (t^ - Zo^)(t - z) 


(D.32) 
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The residuals for expression (D.32) were calculated with equation (D.22). 
They become 


Residuals 



-4I(dz, 


Zo (C + z^)(c - z)| 


+ lim 




-41{?)z, 


(C - z^)(C - z) 


+ lim 
r, z 


-4I(?)z, 


(C - z^){^ z^) 


(D.33) 


To insure single valuedness of the stress functions, values of I(^) 
must be chosen so that they lie on the same branch. Values of 1(c) 
will lie on the same branch for all values for c if 


K?) 

lim = 1 (D.34) 

^->■00 


The two possible values of 1(c) are the complex numbers w^ and -w^. 
Hence, equation (D.34) requires that I(Zq) equals w^ and that I(-Zq) 
equals -w^, or simply that Uz^) = -H-z^). With the proper values 
of 1(c) defined according to equation (D.34), equation (D.33) 


reduces to 
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4S 

Residual = (zl(z.) - z.I(z)) (D.35) 

2 2 _ ^2 0 0 

In order to evaluate equation (D.31) the limits of three 
integrals must be taken, The first two integrals involve limits as 
e ->■ 0 (equations (D.24) and (D.25)). The evaluation of these limits 
is similar for both equations. The limits can be obtained by rearranging 
the terms in the integrand so that they can be expressed by a simple 
binomial series as 

n(n - 1) 

(x + y)" = x'’ + nx^'V + x""^y2...(y2 < x^) (D.36) 

2 ! 


or 


(1 ± x)'^ = 1 + X + x^+ x^ +... (x^ < 1) 


(D.37) 


With the use of the first term of P(t) as shown in (D.32), as an 
example (D.24) can be expressed 


/ A(c) = lim S 
r 2TT 

^ e 0 


.4z„ic - 2aee^®e^®d0 


(ee^® - a + - a - Zo)(ee’® - a - z) 


(D.38) 
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The term Vee^^® - 2aee^® in equation (D.39) can be expressed in terms 
of equation (D.36) as 


- 2a e' 





ee 


i0 




2a 

.10 


ee 


2 2i0 
e e 


4a 


32a' 


(D.39) 


while the term (ee^^ - a + z)”^ can be expressed in terms of equation 
(D.37) as 


1 


1 



a + z 

0 



1 


^0 ■ ^ 


^ 1 



e^e 


2i0 



(D.40) 
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Substituting equations (D.39), (D.40) and similar expressions for the 
remaining terms in the denominator of equation (D.38) into equation 
(D.38) yields an expression as 


1 im 
e 0 


j0 

4z^ 1 e \/^ e ^ 

2TT 



ee 


10 




^0 - ^ 


...j 



(D.41) 

whose limit is zero. Equation (D.25) can be evaluated by similar means. 
The result for P(t) as shown by (0.32) is also zero. 

Equation (D.30) can also be evaluated by a limiting process and 
by the use of equations (0.36) and (0.37) as R ->- «>. For this case, 
the integrand has to be arranged in a manner that makes the expansions 
for equation (0.36) and (0.37) valid for large values of R. For exam- 
ple, for the first term of P(t) as given by (0.37), equation (0.30) 
can be expressed as 
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(D.42) 


By substituting the appropriate expansions for the terms in the integrand 
of equation (D.42), the limit as R increases without bound gives 
the value of equation (D.42) as zero. 

Similarly, equation (D.31) can be evaluated for all ten, is of 
P(t) contained in equation (D.19). All of these terms were combined, 
simplified, and substituted into equation (D.20) to yield 


4>(z,Zq) = n(z,Zg) = S B(z,Zp) + S B(z,Zq) 


(D.43) 


where 


1 r -4z. 


2nz. 


z - z. 


Z + z. 


B(z,z„) = - 


2 iz^-z^" - z^2 (z - Zq)^ (z + Zq)2 


- Z„ 

0 0 


- zz. 


^ -0 1 


2I(z^)I(z) (Zq-z)^ 


I(z) 




Z " - z" 


nl(Zo) 


^0^ - 


(D.44) 


A 
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With the use of equation (D.43), which shows <I>(z) and fl(z) 
to be equal, and equation (D.11), 'I'^z) was expressed as 


M'(z) = $(z) - $(z) - z$'(z) (D.45) 


where o(z) is given by right hand side of equation (D.43). To 
evaluate equation (D.45), the derivative of $(z) was required. Hence, 
the derivative of B(z,z^) was required because 

d4>(z) 

<^‘(z) = = SB'(z,Zq) + SB'(z,¥q) (D.46) 


Using a symbolic manipulation system, MACSYMA (1975), written 
in LISP programming language, the derivative of B(z,z^) was found 
to be 




B'(z,Zp) = -2z 


3z„^z^ + z^z„ - 4z^^ ■) 
0 0 0 0 


(Z^ - Tq^)^ 


(Zo - z)^ (Zq + z)3 




I(Z) |z^ - I (y - 

- 2z“ + a^z^) T 


A 
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(D.47) 

To complete the evaluation of the functions used in equations 
(D.l) through (D.4) for pressure acting on the crack surface, the 
integrals of <i>(z) and 'f(z) given by (p(z) and i1j(z) respectively 
were determined. As shown by equation (D.45), once <})(z) is found 
’p(z) is also determined. To evaluate <))(z) by integrating (D.43) 
the integral of B(z,Zq) was determined (evaluation of the integral 
of B(z,z^) IS the same as for the integral of B(z,z^) except that 
Zg is replaced by z^). 

Many of the terms in B(z,z^) are easy to integrate by using 
standard rules of calculus. However, those terms which contain I(z) 
in the denominator require some rearrangement before the integration 
is attempted. For example, integrals of the type 

dz 

1(z)(Zq - 


Il(z) = / 


(D.48) 
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appear frequently and were integrated as follows. First, a change 
of variables was made by letting 

z = asino, dz = acosedG 

With the restriction that |z| £ a, equation (D.48) became 


acosGdG 


I(z)(Zq - 


I(z)(Zq - asinG) 



At this point care must be taken to choose the correct value of the 
multi-valued function. The correct value was assured by requiring 
equation (D.34) to hold. For example, I(z) can be written as 



but according to equation (D.34) 



= 1 = 1 

Z-voo 2 Z^oo Z 



so then 
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I(z) = - i \/a^ - or \/a^ - z^ = il(z) (D.51) 

Substituting equation (D.51) into equation (D.40) and that result 
into equation (D.49) yields the result 

dz d 0 

/ = 1 / (D.52) 

I(z)(Zq - z) Zq - asinG 

which can be integrated by using standard integral tables. 

Burlington (1973) gives the value of the integral in (D.52) as 


d0 


•1 f-a + ZpSinG +^a^ - z^^cosg' 


Zq - asinG \j 


a^ - zj 


Log 


- sin 0 


(D.53) 


Making the appropriate substitutions of 


sin 0 = z/a 


COS 0 = il(z)/a 
" "o' = 


gives the result of (D.48) as 
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dz 1 f-a^ + zZq - I(Zq)I(z)' 

I(z)(2o - Z) \/2o - 1 

(D.54) 

The development of equation (D.54) was restricted to values of 
]z| < a. However, for values of |z| > a the substitution z = aCSC 0 
can be made and a similar process repeated. The results of this integra- 
tion are Identical to equation (D.54). Consequently, equation (D.54) 
is valid for all values of z (the same can be shown to be true for 
all values of z^). 

After much labor and simplification the integral of B(z,Zq) 
was found to be 


1 


BI(z,z ) * - < Log(z + z.) - Log(z„ - z) + n SLog(z. - z) - Log(z + z 


■ 1 ' 


2(^0 ■ " o ^ 


z" - zj 


Hz) 1 

\Z - Z >+ XI(z,z ) - nXI(z,z )) 

0 J(- J 0 0 

0 _ 


(D.55) 


where 


XI(z,Zq) = Log 


(zz^ - a^ - I(z)I(z^,))(z^ + a)-^ 


(zZq + a^ - I(z)I(z^))(Zj, - z) 
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Therefore, integration of equation (D.43) can be expressed as 

(j)(z. 2 ^) 801 ( 2 , Zq) = Mz,z^)dz => Mz,z^)dz = SBI(z,Zjj) + SBI{z,Zq) 

(D.56) 

By integrating equation (D.45), i|^(z) can be expressed in terms of 
equation (D.56) as 

ij;(z) = ^(z) - z$(z) (D.57) 

Therefore, in summary, the functions required to evaluate 
equations (D.l) through (D.4) for the cracked metal sheet which has a 
stress applied to the crack surface (equal in magnitude but of 
opposite sign to the stress along the crack line for a solid sheet) 
are given by the equations (D.43), (D.45), (D.56), and (D.57). 

Green's Functions for Stress 

As mentioned previously, the solution for the cracked sheet is 
obtained from superposition of the stresses from the two loading 
conditions shown on figure D.la and D.lb. The stress functions 
for these two loading conditions, which were developed in the previous 
sections, were used to obtain the stresses in each loading condition. 

The stresses for the two conditions were then added to form the Green's 
function for stresses for the loading condition shown in figure D.lc. 
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The stress state for four point loads acting on a solid sheet 
shown on figure D.la was found by substituting equation (D.6), 

(D.7), and (D.8) into equations (D.l), (D.2), and (D.3). The result 
is, in terms of the coefficients of the X and Y point loads: 

= ReaHGl - G2 - G3)X + ReaHGlA - i(G2 - G3)}Y 

Oy = ReaHGl + G2 + G3}X + ReaHGlA + i(G2 - G3)}Y (D.58) 

= Img {G1 + G2 + G3}X + Img {GIA + i(G2 - G3)}Y 

xy 

where 


Real 

G1 = 

2it( 1 + n)t^ 




(D.59) 



A 
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1 

G3 = 

2tt( 1 + n)t 
m 




i— • 



(D.62) 

The stress state for the stress applied to the crack surface, 
as shown on figure D.lb, was found by substituting equations (D.43), 
(D.45), and (D.46) into equations (D.l), (D.2), and (D.3). After 
several algebraic manipulations, the result was found to be 


0 ^ = Real{2G5 - G7} X + Real{2G6 - G8}Y 
0 y = Real{2G5 + G7} X + Real{2G6 + G8}Y 
= Img (G7) X + Img (G8) Y 

Ajr 


where 


G5 = 
G6 = 

G7 = 

G8 = 


1 

2tt(1 + 
i 

2tt( 1 + n)t^ 


B(z,Zq) + B(z,Zq) 
B(z,Zq) - B(z,Zq) 


(z - z) 
2tt( 1 + n)t,^ 
id - z) 
2tt(1 + 


B'(z,Zq) + B'(z, 

■ 

B'(z,Zq) - B'(z,Zq) • 

. 



(D.63) 


(D.64) 

(D.65) 


(D.66) 


(D.67) 
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Superimposing the stresses from the two loading conditions by 
adding equations (D.58) and (D.63) and taking the coefficients of the 
X and Y loads for the different stresses yielded six Green's 
functions for stress as 


GS 

GS 


GS^^ = ReaHGl - G2 - G3 - 2G5 + G7} 

(D.68) 

GS = ReaHGlA - i(G2 - G3) - 2G6 + G8} 
xy 

(D.69) 

GS^^ = ReaHGl + G2 + G3 - 2G5 - G7> 
yx 

(D.70) 

GS = ReaHGlA + i(G2 - G3) - 2G6 - G8} 

J 

(D.71) 

(xy)x ’ ^"’9 + G2 + G3 - G7)} 

(D.72) 

(xy)y ~ {GIA + i(G2 - G3 - G8)} 

(D.73) 


where the first index on GS indicates the stress and the second 
indicates the load responsible for it. For example, GS^^ indicates 
the Oy stress at point z due to a unit load applied in the x- 
di recti on at point z^. 

The Green's functions given by equations (D.68) through (D.73) 
were verified with the use of a finite element program developed 
by Y. K. Cheung and I. P. King and documented in Zienkiewicz's book 
(1971). The finite element model used for the test case is shown in 
figure D.3. Because of symmetry only the first quadrant of the cracked 
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Fig. D.3. Finite element mesh used to check Green's functions 
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sheet was modeled; the crack was simulated by freeing the nodes in the 
y-direction along the x-axis from the origin to x = 1.5 inches. 

Point loads of X = 1 and Y = 2 were applied to the model at the 
point Zq = 2.5 + 2.5i, 

The finite element results were compared to the Green's functions 

results. For the comparison, the stresses were calculated along the 

line z = X + 1.251 by equations (D.68) through (D.73) and with finite 

elements. Finite element values of stresses along this line were 

taken as the average of two elements midway between y-coordinates of 

the nodes. On figure (D.4) the dotted line, dashed line, and solid 

line Indicate a , o , and a stresses respectively, obtained 

X y xy 

with the Green's functions while the symijols represent stresses obtained 
from the finite element solution. The comparison was good and verified 
the Green's functions within the accuracy of the finite element model. 

Green's Functions for Displacements 

The displacement field for the solid sheet under four point 
loads as shown on figure D.la was found by substituting the stress 
functions ${z), 4 >(z), and ,|;(t) shown in equation (D.8), (D.9), 
and (D.IO) respectively into equation (D.4). The result is 


2G(u + iv) = SFp(z,Zq) + SFq(z,Zq) 


(D.74) 


where 




Fig. D.4. Verification of Green's functions for stresses in a 
cracked isotropic sheet 
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Fp(z.Zjj) = 2^nReal jLog(z + z^) - Log(z - z 


•'} 


+ 



(0.75) 


The displacement field for stress on the crack as shown in figure 
D.lb was found by substituting the equations for the stress functions 
<()(z), ip(t), and 4>{z) shown in equations (D.56), (0.57), and 
(0.43) into equation (0.4). The result is 


2G(u + iv) = SFjj(z.Zq) + SFq(z.Zq) 


(0.76) 


where 


Fp(z,Zp) = -n BI(z.Zq) + BI(z,Zjj) - (z - z)B(z,Zq) 


(0.77) 


Superimposing the displacement equations for the two preceding 
equations yields the displacement field for the load condition shown 
on figure 0.1c as 


2G(u + IV) = S<Fp(z.Zq) - F 


+ S 


|fp(2 


(0.78) 
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Equation (D.78) was simplified, after lengthy algebraic manipulations. 
The result was found as 


where 


with 


26{u + iv) * SgCz.Zp) + Sg(z.ZQ) 


(D.79) 


n|R( 


g(z,Zp) = n<Re XA(z,Zq) - XC(z 




+ 0.5^ti^XB(z,Zq) + XB(z,Zq)V + XC(z,Zq) 


■"o>] 


- z„z 
0 0 


+ 2 


z^ - z * 

0 J 


^ + (z - z)B(z,Zq) 


XA(z,z.) = - Log <- 


zZq - a* - I(Zq)I(z)' 
zZp + a^ - I(Zq)I(z) 


XB(z,z„) = Log 


zZo - a* - I(Zo)Kz) rz^Zoi 


“o * ■ I(Zo)Xz) I ^0 ■ 


0 0 

XC(z,z„) = 

0 2 2 
Z - 


I(z) 


z - z 


0 T/r 
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1 r -4z, 


2nz, 


B(z.z^) = 




2 /z* - z * z* - z/ (z - zj^ (z 


^^0 ■ 
+ Z )2 

0 J 


z« ■ Z- \ a* - zz„ + zz “ 

0 0 1 0 0 


2I(Zq)I(z) { (z - (z + z^)^ 


I(z) 



n 



■] 


The Green's functions for displacements were obtained from 
equation (D.79) by forming coefficients of the X and Y loads 
for the u and v displacements. The result was 


with 


Real(g(z,ZQ) + g(z,ZQ)) 

(D.80) 

" ‘^o Real (i(g(z,ZQ) - g(z,z^))) 

(D.81) 

'^yx " ^0 

(D.82) 

DGyy = Cq Img (i(g(z,ZQ) - g(z,ZQ))) 

(D.83) 


1 


4Gt tt( 1 + n) 
m 
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where the first index of GD indicates the displacement and the 
second indicates the load responsible for it. 

The Green's functions for displacements given in equations 
(D.80) through (D.83) were verified with the finite element model 
discussed in the previous section and shown on figure D.3. A 
comparison of the displacements calculated from the Green's functions 
and the finite element solution is shown on figure D.5. The comparison 
is made along the line z = x + 2i where the finite element 
displacements are taken from the nodal points of the model. In the 
figure the solid and dotted lines indicate the u and v displacements 
calculated with the Green's functions while the symbols represent 
values obtained from the finite element solution. The agreement is 
good and verifies the Green's functions within the accuracy of the 
finite element model. 




APPENDIX E 


GREEN'S FUNCTION FOR AN UNCRACKED 
ORTHOTROPIC SHEET 

Lekhmtskii (1968) gives the stress and displacements in an 
orthotropic sheet in plane stress in terms of two stress functions, 
$j(Zj) and $ 2 (zJ, as 


Oy = 2 Real^$j'(Zi) + 


a = -2Real 


y + u. 


= 2Real jsj + s/<I )2 (z^ij 

|^$i'(Zi) + ^>2'(z2)j 
|si$i'(Zi) + S2$2' (Z2)j 

u = 2 Real|pj(i,j(Zj) + P2$2'(^2)j " 

V = 2Real|qj$^(zJ + q2«I>2(z2)| + w^x + v^ 

where Sj and Sj are the roots of the equation. 


CE E 

X X 

s'. + + — = 0 

G E 

t xy _ y 


(E.l) 

(E.2) 

(E.3) 

(E.4) 

(E.5) 


(E.6) 
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where E , E , 6 , and v are the elastic moduli and Poisson's 
A y xy xy 

ratio. Leknitskii (1968) proved that the roots, s^ and s^, 
could not be purely real for real materials but are purely imaginary 
or complex. The complex roots occur in conjugate pairs and s^ and 
s^ are distinct roots, i.e. not complex conjugates. The following 
were defined 


:^ = X + Sjy = X + s^y 




yx 


1 1 

P. = 7 - V) P, = 7 - \y) 

X "^x 


(E.7) 


= 


^'y 


(1 - V s M q = (1 - V s M 

' yx 1 ' ^2 . n ' yx 2 ^ 


s,E. 


yx 2 


(E.8) 


Note that w^, u^, and in equation (E.4) and (E.5) are rigid 

body rotations and translations respectively and 



dfj(zJ 


dz. 




d^jCZj) 


dz. 


The two required stress functions for a point load acting on a 
solid, orthotropic sheet were given by Lekhnitskii (1968) as 


4j(z,) = A^Logz, 


= Bj.Logz^ 


(E.9) 
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where A and B satisfy the equations (for material axis concurrent 
^ c 

with the axis of orthotropicity) 






Y 


2irt^i 

c 


s + s - s - s = 
1 c 2 c I c 2 c 


ZTTt^i 


s ^A_ + s ^B_ - s ^A_ - s ^B_ 
1 c 2 c 1 c 2 c 


-Ay'' 


ZlTt^l 

C 



Vy/X 


2-rrt_i 

C 


(E.IO) 


With the use of the preceding equations, stress functions for 
point loads acting on a unidirectional boron/epoxy composite were 
developed. For this material the material constants are 


= 0.27 X 10^ psi E„ = 3.0 x 10^ psi 
X y 


= 0.7 X 10° psi = 0.019 

Ajr Ajf 


With these material constants, the roots of equation (E.6) were found 
to be purely imaginary as 
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s = ±0.1531 s = ±1.951 

1,3 2,4 

For purely Imaginary roots equations, and were found from 
equations (E.IO) as 


where 


= Cn X + 1 C12 Y 


= C21 X + 1 C22 Y 


(E.ll) 

(E.12) 


Cll = 


s (s + 1) 

2 1 xy 

Airt (s 2 - s 

c 2 1 


C12 = 




4iTt (s 2 - s 
c 2 1 


C21 = 


4irtc(s^^ - s^*) 


C22 = - 


s + V y 

2 xy 


4iTt (s ^ - s 

^ Z I 


Using equations (E.ll) and (E.12) and translating the origin 
so that the singularity occurs at the point z^, equations (E.9) 
became 


4>i(Zi,Wj) 


(Cll X+ 1 C12 Y)Log(z - w ) 

1 1 


$ (z ,w ) = (C21 X+ 1 C22 Y)Log(z - w ) 
2 2 2 2 2 


(E.13) 
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where w = x„ + s w = x„ + s 

1 0 1-^0 2 0 2 0 

Equations (E.13) were used to construct a solution for four 
point loads acting on a solid orthotropic sheet as shown in figure E.l. 
The result was 

$ (z ,w ) = Cll G8(z ,w )X + i C12 G9(z ,w )Y 
111 11 11 

(E.14) 

3> (z ,w ) =C21 G8(z ,w )X + i C22 G9(z ,w )Y 
2 2 2 2 2 2 2 


where 


G8(z,w) = Log(z - w) - Log(z + w) - Log(z + w) + Log(z - w) 

(E.15) 


G9(z,w) = Log(z - w) + Log(z + w) - Log(z + w) - Log(z - w) 

(E.16) 

The derivatives of equations (E.14) were found to be 

$ '(z ,w ) = Cll G8'(z ,w )X + i C12 G9'(z ,w )Y 
111 11 22 

(E.17) 

$ '(z ,w ) = C21 G8'(z ,w )X + i C22 G9'(z ,w )Y 
2 2 2 1 1 2 2 


where 




Fig. E.l. Location of point loads on orthotropic sheet 
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G8'(z,w) 


2w 2w 



(E.18) 


2w 2w 

G9'(z,w) = 3 - 

z2 _ Z^ - 


(E.19) 


Equations (E.13) and (E.17) can be used in equations (E.l) 
through (E.5) to describe the stresses and displacements in the 
orthotropic sheet. 


Green's Functions for Stresses 

The Green's functions for stresses were developed by 
substituting equations (E.17) into equations (E.l), (E.2), and (E.3) 
to determine the stress state as 


’’x = «Sxx* " 

(E.20) 

"y = «yxX * 

(E.21) 

“xy = H5(xy)xX * «S(xy)y'' 

(E.22) 


where the Green's functions are given by 
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HS 

HS 




= 2Rea1<s ^ C11 68' (z , w ) + s ^ C21 68' (z , 

XX I 1 112 2 


HS^y = 2Rea1 




«,)j 


2 C12 69' (z ,w ) + s C2269' (z ,w ) 
112 2 2 


' 1 ! 


= 2Real< C11 68' (z ,w ) + C21 68' (z ,w ) 
yx I 11 2 2 


yy 


(xy)x 


(xy)y 




i 


, 


= 2ReaN i { C12 69' (z ,w ) + C22 69' (z , 


1 1 


= -2Real<s Cll 68' (z ,w ) + s C21 68' (z 


I 


'1 

’".'I 


-2Real 




1 1 2 


C12 69' (z ,w ) + s C22 69' (z 
112 2 




(E.23) 

(E.24) 

(E.25) 

(E.26) 

(E.27) 

(E.28) 


To verify these 6reen's functions, stresses computed with equa- 
tions (D.20) through (D.22) were compared to finite element results. 
The model used for the finite element solution is identical to that 
shown in Appendix D on figure D.3 except that no nodes were freed 
along the x-axis to simulate a crack. The finite element program used 
is documented by Zienkiewicz (1971) as mentioned in the previous 
appendix. However, the program as given by Zienkiewicz does not have 
orthotropic capability. Therefore, it was modified by introducing 
the orthotropic stiffness matrix for plane stress given by 
Zienkiewicz (1971) as 
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1 0 
0 M(1 - nvy^2) 

in place of the isotropic 

The comparisons between the stresses along the line z = x + 

1.25i are shown on figure E.2 on which the dotted line, dashed line, 

and solid line indicate the o , a , and a stresses respectively 

A y xy 

from the Green's functions. The symbols on the figure indicate the 
same stresses obtained from the finite element solution. The 
comparison was good and verified, within the accuracy of the finite 
element model, the accuracy of the Green's functions for stresses. 

Green's Functions for Displacements 

The Green's functions for displacements were developed by 
substituting equations (E.14) into equations (E.4) and (E.5). Because 
of the symmetry of the loads the rigid body rotation, w^ , and the 
rigid body translations, u^ and v^, are zero in equations (E.4) and 
(E.5). The results are 

u = HD^^X + HD^yY (E.29) 


[D] 


1 - n\) 


nv 


yx 


yx 


Ex Sx 

where n = — and m = 

Ey E^ 


stiffness matrix in the program. 


V = HD X + HD Y 
yx yy 


(E.30) 



stress, psi 



Fig. E.2. Verification of Green's functions for stresses in 
an uncracked orthotropic sheet 
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where the Green's functions are given by 




= 2Real (p Cll G8(z ,w ) + p C21 G8(z ,w ) 

XX J 1 112 22 


{'h 


= 2ReaW i {p C12 G9(z ,w ) + p C22 G9(z , 

xy ill 112 2 


HD,„ = 2Real< q Cll G8(z ,w ) + q C21 G8(z ,w ) 

yx J 1 112 22 


HDyy = ?Real 




C21 G9(z ,w ) + q C22 G9(z , 
112 2 


'] 

v]! 

■41 


(E.31) 


(E.32) 


(E.33) 


(E.34) 


With the use of the finite element model shown in figure D.3, 
equations (E.31) through (E.34) were verified by comparing the displace- 
ments along the line z = x 2i calculated with equation (E.29) 
and (E.30) shown on dotted and dashed line, respectively, on figure 
E.3 with the displacements calculated with the finite element program. 
The comparison is good and within the accuracy of the finite element 
model verifies the Green's functions for displacements in the 
uncracked orthotropic sheet. 




APPENDIX F 


COMPUTER PROGRAM TO PREDICT CRACK 
AND DEBOND GROWTH 

The analysis developed in Chapters IV and V to predict crack 
growth in reinforced systems requires the use of a digital computer 
to perform numerical integration and solve large systems of simultaneous 
equations. The analysis was programmed in FORTRAN IV for use on the 
NASA Langley Research Center (LRC) CDC 6600 computing system. However, 
no special system routines were used in the program, and the program 
should be usable on any computing system that uses a FORTRAN IV 
compiler. With the exception of the Gaussian elimination subroutine 
SIMQ, the program is all original code. 

For economical reasons the user should be aware of the central 
processing time (CP) and central memory (core) required to execute 
the program. The CP time is primarily a function of the numerical 
integration of the Green's functions. Several integrations are 
required for each shear element, i.e. elements of region B shown on 
figure 8. Consequently, the CP time requirement is a function of the 
number of shear elements and can vary from a few to several thousand 
seconds depending upon the number of elements. The core requirements 
are also a function of the number of shear elements. Each element 
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contains two unknown shear stresses, and The number of 

rows and columns of the square, fully populated coefficient matrix 
used to solve for these stresses is two times the number of elements. 
In the program the dimensions are set for fifty shear elements 
(100 unknowns), but the program can accomodate more elements if the 
array Z9(10,000) in the program is enlarged. For fifty elements the 
core requirement is 137K octal. 

Data are read in the program via a NAMELIST with the following 
definitions : 

E3 modulus of the cracked sheet 
T2 thickness of the cracked sheet 
VI Poisson's ratio of the cracked sheet 
E4 modulus of the composite sheet in the y-direction 
(loading-axis) 

E5 modulus of the composite sheet in the x-direction 
V4 Poisson's ratio of the composite sheet for a load applied 
in the x-direction 

GC shear modulus of the cracked sheet 
T4 thickness of the reinforcement sheet 

TAD thickness of the adhesive layer 

GAD initial effective shear modulus of the adhesive 
GAD2 secondary effective shear modulus of the adhesive 
SYIELD yield stress of the bulk adhesive in uniaxial tension 
A1 initial crack length in metal sheet 

F initial aspect ratio of the debond ellipse 
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NC number of columns of interlaminar sheat elements 

NR number of rows of interlaminar shear elements 

FCYC final number of applied load cycles 
S remote applied stress applied to the reinforced system 
in the y-direction 

As an example a sample run was made for one iteration. Figure 
F.la shows the model prior to the iteration and figure F.lb shows 
it after the iteration. The sample NAMELIST input for the run was 
as follows: 



s = 25,500 psi 


I I I I I I 



s = 25,500 psi 


y M I I I I 

I 



00 

CO 


Fig. F.l. Model for sample run 
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SAMPLE INPUT 

$INPU E3=1.0E+07,T2=.156,VI=.30,E4=3.0E+07,E5=3.0E+06,V4=.02,GC=7.0E+05, 
TAD=.004,GAD=65000. ,GAD2=36000. ,SYIELD=4200. ,A1=1 .00,F=.001 ,NC=5,NR=3, 
FCYC=10000. ,S=25500. ,T4=.0208, 

$ 


SAMPLE OUTPUT 


MODULUS OF THE CRACKED SHEET .107E+08 
THICKNESS OF THE CRACKED SHEET .156E+00 
POISSONS RATIO FOR THE CRACKED SHEET . 300E+00 
COMP MODULUS PARALLEL TO LOAD .300E+08 
COMP MODULUS TRANSV TO LOAD .300E+07 
SHEAR MODULUS OF COMPOSITE .700E+06 
COMP POISSONS RATIO TRANSV TO LOAD AXIS .200E-01 
THICKNESS OF THE REINFORCEMENT SHEET .280E-01 
MINOR AXIS HEIGHT IN PER CENT CRACK LENGTH .lOOE-02 
REMOTE APPLIED STRESS .255E+05 
NUMBER OF COLUMNS FOR BOUNDARY POINTS 5 
NUMBER OF ROWS FOR BOUNDARY POINTS 3 
FINAL NUMBER OF APPLIED LOAD CYCLES .lOOE+05 
INITIAL CRACK LENGTH .100E+01 
ADHESIVE THICKNESS .400E-02 
SHEAR MODULUS OF THE ADHESIVE .650E+05 
SECONDARY SHEAR MODULUS AFTER YIELDING .360E+05 
YIELD STRESS OF MODULUS IN UNIAXIAL TENSION .420E+04 


METAL 

COMPOSITE 


SIGM-X= 

SIGM-X= 


.269E+03 
- . 1 50E+04 


SIGM-Y= 

SIGM-Y= 


.200E+05 

.559E+05 


X AND Y DIMENSIONS OF INTEGRATION AREA ARE . 990E+00 


.746E+00 
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NODE NUMBER 

X-COR 

Y-COR 

1 

.990E-01 

.125E+00 

2 

. 297E+00 

.121E+00 

3 

.495E+00 

.113E+00 

4 

•693E+00 

•988E-01 

5 

.891E+00 

.765E-01 

6 

.990E-01 

.373E+00 

7 

. 297E+00 

.365E+00 

8 

.495E+00 

.349E+00 

9 

. 693E+00 

•323E+00 

10 

.891E+00 

. 285E+00 

11 

.990E-01 

.fa22E+00 

12 

.297E+00 

•613E+00 

13 

.495E+00 

.593E+00 

14 

.693E+00 

.563E+00 

15 

.891E+00 

.521E+00 

14 ELEMENTS HAVE YIELDED AT 

•189E+05 



YIELDED ELEMENTS 





3 2 

1 

4 

5 

9 

8 7 

10 

6 




THE 

YIELD MACROSCOPIC STRESS 

IS .716E+04 

ELASTIC 

NODE 

DEBOND HEIGHT 

X-COEFFICIENT 

Y-COEFFICIENT 

1 

.125 

.300E+01 

.849E+04 

2 

.121 

.279E+02 

.859E+04 

3 

.113 

.lOOE+03 

.864E+04 

4 

.099 

.327E+03 

.847E+04 

5 

.076 

.lllE+04 

.667E+04 

6 

.373 

.346E+02 

.256E+04 

7 

.365 

.130E+03 

.266E+04 

8 

.349 

.328E+03 

.276E+04 

9 

.323 

.806E+03 

.276E+04 

10 

.285 

.191E+04 

.220E+04 

11 

.622 

.916E+02 

.117E+04 

12 

.613 

.293E+03 

.119E+04 

13 

.593 

-560E+03 

.117E+04 

14 

.563 

•983E+03 

.108E+04 

15 

.521 

.180E+04 

.807E+03 
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NODE DEBOND HEIGHT 


WITH PLASTICITY 

X-COEFFICIENT Y-COEFFICIENT 


1 

.125 

.207E+01 

.708E+04 

2 

.212 

.213E+02 

.709E+04 

3 

.113 

.795E+02 

.704E+04 

4 

.099 

.245E+03 

.679E+04 

5 

.076 

.805E+03 

. 536E+04 

6 

.373 

.438E+02 

.312E+04 

7 

.365 

.156E+03 

.321E+04 

8 

.349 

.371E+03 

.325E+04 

9 

.323 

.864E+03 

.315E+04 

10 

.285 

.191E+04 

.244E+04 

11 

.622 

.105E+03 

.165E+04 

12 

.613 

.333E+03 

.165E+04 

13 

.593 

.635E+03 

.161E+04 

14 

.563 

.110E+04 

.146E+04 

15 

.521 

.189E+04 

.112E+04 

K-BOND 

K-UNSTIFFEND 

K-STIFFENED 

K- FACTOR 

-.108E+05 

.200E+05 

.919E+04 

.459E+00 


STRESS AND STRAIN 

IN THE METAL 



NODE 

SIG-1 

SIG-2 

SIG-12 

1 

6 

11 

.681E+04 
. 340E+04 
.137E+04 

.580E+04 

.121E+05 

.154E+05 

-.472E+02 

-.103E+03 

-.146E+03 


NODE 

EPS-1 

EPS-2 

EPS-12 

1 

-.799E-03 

.733E-03 

-.115E-04 

6 

-.657E-03 

.123E-02 

-.250E-04 

11 

.559E-03 

.148E-02 

-.354E-04 
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STRESS AND STRAIN IN THE COMPOSITE 


NODE 

SIG-1 


SIG-2 

SIG-12 

1 

.838E+04 


.105E+06 

.742E+03 

6 

.897E+04 


.728E 05 

.252E+03 

11 

.758E+04 


. 568E+05 

-.301E+03 

NODE 

EPS-1 


EPS-2 

EPS-12 

1 

.272E-02 


.350E-02 

.674E-04 

6 

.294E-02 


.242E-02 

.360E-04 

11 

.249E-02 


.189E-02 

-.430E-04 

ENERGY RELEASE ON 

MINOR AXIS 

IS 

.597E+01 


APPLIED CYCLES= 

. 990E+01 




DA/DN BEFORE CYCLE 

INCREMENT 

IS . 

181E-04 


DB/DN BEFORE CYCLE 

INCREMENT 

IS . 

202E-01 


INCREMENT CYCL = . 

990E+01 

CRACK 

,LENGTH= .lOOE+01 



MAX DEBOND HEIGHT= .201E+00 
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HIOGRAM COLLS( INPUT .OUTPUT .TAPE4 ) 

C 

C THIS PROGRAM CALCULATES CRACK AND DEBOND GROWTH IN A 
C REINEORCED SYSTEM COMPOSED OF A CRACKED METAL 
C SHEET REINFORCED VHiTH A UNIDIRECTIONAL BORON/EPOXY 
C SHEET. 

C 

COMMON /ROT/Z9( 100 00),D( 100 ),NALF 

OOMMON/BOND2/NE,NL,NT,XD( 100),YC( 100),XA( 100 ),YA( 100 ), 
INOP(IOO) 

DIMENSION FTNC(50),PXNCY(50),FA1(50),FB1(50 ) 

COMMON /ADHES/TAD , GAD 

COMMON /TO P/E3, T 2. VI , SMX , a^IY , G, CONS . Q, A1 , SCON 
DIMENSION DR( 100 ) ,CPR( 100 ) ,DPR( 100 ) 

COMMOJI/XLIMIT/DC( 2,51 ),DB(6, 51 ),EP 
DIMENSION CM(3,5),SIG(5),STRAIN(3),DD(100 ) 

DIMENSION SMbi3),a*IM(3,3),CMC( 5, 3) ,STRESSM( 3), STRESSC 
1(3) 

DIMENSION TSIGM(3),TSIGMT(3),31M(3,3),31C(3,3) 
DIMENSION TOTG( 50 ) , SIGMT( 3 ) , TSTRA3N( 3 ) , STRESS ( 2 , 3 , 50 ) 
1,STRANN(2,3,50) 

COMMON /BO T/E4 , T4 , V4 , sex , SCY , GE , OONB , Q 1 , GC , E5 

COMMON/CT 0 L/TO L , NC , NR , TX , TY , NBC ( 100 ) , IBC 

COMMON /BO ND/F , P 1 , P2 , XKF , XKUN S , XKS T IF 

COMMON /DEC AY/EE ,ALPHA1 , ALPHA2 

00MPI2X Z,CI,XK 

EXTERNAL XK 

CI=CMPLX(0. ,1.0) 

NAMELIST /INPU/E3, T2, F , 3, TEST , T4, E4, V4, VI , E5, GC , TAD , GAD 
1,NR,NC ,FCYG, 

1A1,GAD2,SYIELD 

C 

C E3-M0DULUS OF THE CRACKED SHEET 

C T2-miCKNESS OF THE CRACKED SHEET 

C VI- IS THE POISSONS RATIO FOR THE CRACKED SHEET 
C GC -SHEAR MODULUS OP THE REINPORCIMENT SHEET 
C E4- MODULUS OF THE REINFORCEMENT SHEET IN THE LOADING 

C DIRECTION 

C E5- MODULUS OF THE REDTP ORCEMJINT SHEET TRANSVERSE TO 

C THE LOADING AXIS 

C V4 IS THE POISSaiS RATIO FOR THE BOTTOM SHEET FOR A 
C TRANSVERSE LOAD 

C T4-'IHICKNESS OF THE REDTPORCEI®NT SHEET 

C F-MTNOR AXIS HEIGHT IN PERCENT CRACK LENGTH 
C S -REMOTE APPLIED STRESS 

C NR -NUMBER OP ROWS OF BOUNDARY lOINTS 

C NC -NUMBER OF COLUMNS OF BOUNDARY POINTS 

C ECYC-FINAL NUMBER OP LOAD CYCLES 

C A1-DTITIAL CPJLCK LENGTH 

C 

K=0 



o o o 
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1000 CONTINUE 
READ INPU 
MX=NR»NC 
M=MX 
N=MX 
NALF=MX 

OUTPUT IHE INPUT 

PRINT 1,E3 
PRINT 2,T2 
PRINT 20, V1 
HUNT 18, E4 
PRINT 26, E5 
PRINT 27, GO 
HUNT 21, V4 
PRINT 19, T4 
PRINT 3,F 
PRINT 7,S 
PRINT 30, NC 
PRINT 31, NR 
PRINT 32,FCYC 
HUNT 34, A1 
HUNT 24,T/ID 
HUNT 25, GAD 
PRINT 5,GAD2 
PRINT 6,SYIEID 
P1=2. 

P2=2. 

T0L=.001 
FF = F 
N=2»^N 
M= 2*K 
NSQ=N*N 

IF(N.GT.100)500, 501 

500 HUNT 502, M,N 

502 FCRMAT( //•*< ERRORS — M CR N EXCEEDS DIMENSIONS BOUNDS 
1K=*I10*N=*I10//) 

GO TO 1001 

501 CONTINUE 
Q=(3.-V1)/(1.+V1) 

Q1=(3.-V4)/(1.+V4) 

G=E3/(2.*h .+V1)) 

C0NS=1./(12.566*T2 *(1-+Q)*G) 

SC0N=1 ./(6.2832*( 1-»0 )*T2) 

GB=E4/(2.*( 1.+V4)) 

C0NB=1./(12.56 6*T4 *(14Q1)*GB) 

V5=V4*E5/E4 

EE=2. 7182818 

ALPHA 1 =GAD»( 1 . / (T 2*E 3 > 1 . / (T 4 *E 4 ) ) /TAD 
ALPHA 1=SQRT( AUHA 1 ) 



oooooo o o o o o o ooo o o o o o 
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ALPHA2=GAD2^( 1 . / (T 2*E 3 )+ 1 . / (T 4*E 4 )) /TAD 
ALPHA 2=SQRT(ALtm2) 

TY=-A10G( .05 )/ALPHA2 

CALCULATE THE REMOTE STRESSES 

ZERO ALL COMPLIANCE MATRICES 

DO 50 1=1,3 
DO 50 J=1,3 
a^i(l,J) = 0. 

CMM(I,J) = 0. 

a^c(i,j)=o. 

50 CON T HUE 

CALCULATE COrPLIAIICE FOR METAL SHEET 

Cia^(l,l)=E3/(l.-V1^-2) 
a-IM(1,2)=V1KlMM(l, 1 ) 

CI/IM(2,l)=CIFi(l,2) 

CMM(2,2)=CM(1,1) 

CMM(3,3) = G 

GENERATE UIE STIFF MATRIX FOR THE METAL SHEET 

CALL INVER ( Q-IM , , DD ) 

CALCULATE Ca-IPLIANCE FOR CaiPOSITE 3IEET 

:{N=E5/E4 
.<!-I=GC/E5 

CC=E4/(1.-XN*V4*^ 2) 

CT^C(1 , i}=:oi*cc 
ci-ic(i, 2)=:ai*v4*cc 

CMC(2, 1)=CI-E( 1 ,2) 

CMC (2 2) = CC 

CT'IC ( 3 ! 5 ) = XI > ( 1 . -XN-^ V4->^ 2 )* CC 

GENERATE OHE STIFF MATRIX FOR THE COMPOSIDE SHEET 

CALL INVER(CMC,S'IC ,I!D) 

GENERATE FiACRO STIFFNESS mTRIX 

DO 51 1=1,3 
DO 51 J= 1 , 5 

Cl': ( I , J) = CI'M( I , J) * T2+CMC (I , J) =< T4 

51 OONTHUE 
C 

CALL INVER(CM,ai,I!D) 
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3IG(1)=0.0 

3IG(2)=S*(T2+T4) 

SIG(3)=0.0 

CALL MULT ( SM, SI G, STRAIN) 

400 K)RMAT( 3E11 .3) 

C 

C CAICULATE STRESSES IN EACH LAYER 
C 

CALL MULT (CMM. STRAIN, SIRE SSM) 

SMX=STRESSM(1) 

3‘'IY=STRESSM(2) 

CALL MUHr(CMC.STRAIN,STRESoC) 

SCX=STRISSC(1) 

SCY=STRSSSC(2) 

IRINT 300 ,SI-iX,^IY 

300 FORMATC /i* META.L SIGM-X=*E 10 . 3" SIGM-Y=*E10. 3 ) 

PRINT 30 1, sex, SC Y 

301 P0RMAT( A COMIOSITS SIGM-X=*E10.3- SIGM-Y=*E10.3/) 

c 

C SET Ur LOOIS FOR RIGHT EAW VECTOR AND COEFFICIENT MATRIX 
C 

T1=T2 
CALL CPAR 
C 

C LOOP ON IIlCREKEtTTS OP LOAD CYCLES 
C 

INC =0.0 
K0UNT=0 
1004 CONTINUE 

IF(i:0UNT.IIE.0)]:RIIIT 551 ,KOUNT 

551 FORI''AT(//:OX* XXXXXXXIffiXX UICREI-ENT FOR LOAD CYCLES^ 
1->‘IS->I5 XXXXXXXXX^/) 

TX=.59^A1 
PRINT 16,TK,TY 

16 P0RI:AT(A X ANDY DTMINSI®S OF INTEGRATION AREA ARE ^ 
12E11 .3/) 

E0UT]T=K0UNT+1 

MSAVE=i: 

CALL GRED 
CALL PORI.(N,M) 

C 

C MAiOS RIGHT HAI7D SIDE A UNIT VECTOR FOR USE \/ITH PLASTIC 
C AI'IALYSIS 
C 

DO 1401 1=1, N 

C SAVE UNIT RIGHT II/JID VECTOR 
IH( I)=D(I)/S 
D(I)=D(I)/S 
1401 CONTINUE 
REWIND 4 

1130 F0RINiT(l1P10.3) 
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C 

C STORE 00SPITCIEI:T MA.TRIX OH TAPE 

c 

;^:iITE(4, 1130 )(Z 9(1) ,1=1, 1'JSQ) 

rev; I in) 4 

CALL SIHQ(Z9,D,K,niD) 

c 

C RETRIVE OOEFEICISIIT FiATRIX 
C 

RSAD(4,1 130)(Z9(I) ,I=1,NSQ) 

C 

C SAVE ELASTIC SOLUTia: 

c 

DO 1402 1=1,11 
DD(I)=D(I)-S 
1402 CONTIIRJE 
C 

C CALL THE PLASTIC SUBROUTINE TD FIND INTERLAIIINAR STRESSES 
C AITER YIELDING 

C 

CALL PIA3TIC(SYIELD,GAD2,S,ni) 

NALItN/2 
IRINT 17 

17 K)RKAT(/70X’<V/ITH i'LASTICITY*293P^ELASTIC^ ) 

PRINT 11 
DO 104 I=1,NA.LF 
Ii: = I+lTAIF 
X=DC( 1,1) 

Y=DC(2,I) 

PRINT 12,X,I,Y,D(I),D(IK),IE)(I) ,DD(I+NALF) 

104 CONTINUE 
C 

:<?:tot=o. 

DO 105 J=1,NA.],F 
j';;=DC(i,j) 

Y=DC(2,J) 

Z=X+CI^Y 

CALL XD-ITG(Z,J,A1,S11 , S12 , S21 , S22 ,:<K ,3 ) 

XK1=S11+S12 

aC2=-(S21+S22) 

XKT0T=XKT0T+XK1 

105 OC'NTINUE 
C 

C IRIIIT THE STRESS INTENSITIES 
C 

:(KUN3=SIir-A1>‘->' .5 
.<KS Tn=’= XNUNS + XKTO T 
.KF=XNSTIF/XKUNS 
IRi:iT 14 

i RUT? 15 , X'lTOT ,XICUK3,XKSTIF,XKF 
100 CONTINUE 



198 


M=MSAVS 

C 

C CALCUIATS STRAINS AKD STRESSES AT COORDINATE ROINTS IN 
C THE ADHBREtIDS 
C 

DO 1220 I=1,r>K,IE 
X=DC(1,I) 

Y=DC(2,I) 

Z=X+CI^Y 

DO 12 26 KTYVD=1,2 

CALL VE RI( Z , S(X , SSY , S(XY ,SXX,Sn, SYKY ,MX , KTYIE ) 

C 

C FIRST THRIVE TERI-iS ARE FROM X-L6ADS THE LATTER THREE FROK 
C Y-DOADS 
C 

TSIGH(1)=3XX+SYX 
TSIGri(2)=SXY+SYY 
TS I GH( 3 )= S X S Yr/ 

C 

C CAICUIATE THE STRESSES IN THE COMPOSITE SHEET 
C 

CALL RESIGM(Z,S11 ,S22 ,S12 ,ErYPE,STRESSI!,STRESSC) 

C 

C SUPERIMPOSE STRESSES 
C 

TSI GMT ( 1 )= T SIGM (1 )+S 1 1 
TSia-IT(2)=TSIGM( 2)+S22 
TS IGMT ( 3 ) = T SIGN ( 3 )+ S 12 
C 

C CAI.CULATE TIK CORRESPONDING STRAINS 
C 

IF(i:TY-i^E.E5 .1 )1221 , 1222 
12 21 CALL I-ULT(SIE:,TSIGnT,TSTRAIN) 

GO TO 1224 

1222 GALL mLT ( 3^C ,T3IGMT,TSTRAD: ) 

12 24 CONTDIUE 

c 

C STORE ALL VAIUES 
C 

DO 1225 E:=1,5 

3TR ESS ( I:T Y1>E , DC , I) = T SIGMT ( IX ) 

STR Aini( HTYPE , D1 , 1) =T STIUIN(IK ) 

1225 CCNTIIRJE 
12 26 CCKTEIUS 
1220 CONTINUE 
PRINT 1227 

12 27 FOPJIATC/* stress AlH) STRAIN IN THE METAL ^/) 

PRINT 1228 

1228 F0RI'1AT( 2X>^rOIE- 5:'P< SIG-1 *8X>^SIG-2*8X- SIG-12 *8X^£PS-1 *8X 
1*EPS-2 CX^ EPS-12*) 

PRINT 1229 , (I, ( STRESSd ,K,I) ,K=1,3) , (STRANNd ,K,I) , 
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1K=1,5),I=1,MX,NC) 

1229 FOR1'1AT(I5,3XS10.3,3XE10.3,3XE10.3,4XE10.3,3XE10.5,3XE1 

10.3) 

PRINT 1230 

1230 FORMATC/* stress and SmAIN IN THE COMPOSITE*/) 

PRINT 1228 

PRINT 1229,(1 .(STRESS(2,K,I) ,K=1,5) , (STRANN(2,K,I) ,Ifc1 

1.3) ,I=1,MX,NC) 

SUM=0. 

DX=TX/NC 

DO 1232 I=1,MX,NG 
DY=TY/NR 

SUM=SUH+DY* ( STRANN (2,2,1) -STRANN ( 1 , 2 , 1) ) *D(MX+I ) 

1232 CONTIIUE 

PRim i25i,sroi 

1251 F0RI'^AT(/* energy release on minor axis IS*E11.3/) 

CALL XINC(SUM,MTCY,D^S,I5'S) 

CPR(KOUNT)=DAS 

DPR(K0UNT)=DES 

IF(XKSTIF.(E.56000.)G0 TO 1005 
TNC=TNC+YNCY 

IF(ONC.C2).FCyC)GO TO 1003 

B1=F*A1 

HUNT 42, BTC 

42 K)RMAT( /* APPLIED CYCLES=*E10 . 3) 

HUNT 555,DAS,DPS 

555 FORMdT(* m/IN BEFORE CYCLE INCREMENT IS*E11 .3/ 

1* DB/M BEFDRE CYCLE INCREMENT IS*E 11 .3/) 

PRIITT 40,HTCy,A1,B1,P1,P2 

1* MAXIMUM DEBOND HEIGHT=*E10.3* P1=^»E10.3* P2=*E10.3/) 
40 FORMAT( 5X* INCREMENT CYCL=*E10 . 3* CRACK ,LENGTH=*E 10 . 3 
FTNC(KOHTT)=TNC 
FXNCY(KOUNT)=XNCY 
FA1(K0UNT)=;A1 
FB1(K0UNT)=B1 
rP(A1.GT.2)G0 TO 1001 
GO TO 1004 

1 F0RI-I/^T(20™0DULUS OF THE CRACKED SHEET 

1 *,E12.3) 

2 FORMAT (20»*THICKNESS OF THE CRACKED SHEET 

1 *,E12.3) 

3 F0FI'IAT(20:^MIN0R AXIS HEIGHT IN PERCENT CRACK LENGTH 

1 *,E12.3) 

4 FDRMAT( 20 INCREMENTS TO MXIUMUM CRACK LENGTH 

1 *, 112 ) 

7 FORMAT ( 20 REMOTE APPLIED STRESS 

1 »,E12.3) 

8 F0RMAT(20X^INITIAL CRACK LENGTH BEFORE CYCLING 

1 *E12.3) 

9 FORT'IAT(/25X^INCREMENT ^<I5,4X*CRACK LENGTH *F10.3/) 

10 F0RMAT(5E12.3) 
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11 POm^!AT(203?^LOCATION NODE DEBOND HEIGHT*6X 

1*X-COEPFICIENa>»^8X^Y-OOEPPICIENT*73P^X-COEFPICIENT* 
27X»Y-GOEPPIGIENT*/) 

12 PORMAT(22XP5.3,9XI2,10XP5.3,4(10XE10.3)) 

14 P0RI>1AT( /233P^K-B0ND*5iP^K-UNSTIPPEND*52*K-aPIPPENED*5X 
1*K-PACTOR*/) 

15 PORPIAT(20X,3(E10.3,5X),E10.3/) 

5 P0RMAT( 20 30^ SECONDARY SffiAR MODULUS AITER YIELDING 

1 *E12.3) 

6 P0RMAT( 20 30^ YIELD STRESS OP MODULUS IN UNIAXIAL TENSION 

1 *E12.3) 

18 P0RMAT( 20 30^0 OMP MODULUS PARALLEL TO LOAD 

1 *,E12.3) 

19 P0RMAT(2OX*TIIICKNESS OP THE REINPORCEMENT SHEET 

1 *,E12.3) 

20 P0RHAT( 20 3P^ POISSONS RATIO POR THE CRACKED SHEET 

1 *,E12.3) 

21 P0RMAT(2OX*COI«IP POISSONS RATIO TRANSV TO LOAD AXIS 

1 *,E12.3) 

22 FORMAT( 20X»NUMBER OF ELLIPSES IN GRID 
1 * 112 ) 

23 P0RMAT( 20X*NUHBER OP HNES IN GRID 

1 *1 12 ) 

24 FORM AT( 20 3P^ ADHESIVE THICKNESS 

1 *,E12.3) 

25 FORM AT( 20 3P^ SHEAR MODULUS OP THE ADHESIVE 

1 *,E12.3) 

26 PORMAT(203P^COMP MODULUS TRANSV TO LOAD 

1 *,E12.3) 

27 FORM AT( 20 3^ SHEAR MODULUS OF COMPOSITE 

1 *,E12.3) 

30 FORMAT( 20X*NUMBER OF COLUI-DJS POR BOUNDARY POINTS 
1 *, 112 ) 

31 FORMAT( 203P<NUM:ffi:R OP ROV/S FOR BOUNDARY POINTS 

1 *, 112 ) 

33 FO RMA T ( 20 3C>^ FINAL IHJMBER OP APPLIED LOAD CYCLES 

1 *,E12.3) 

34 FORMAT( 20 3P^ INITIAL CRACK LENGTH 

1 *,E12.3) 

1005 PRINT 38,THC 

38 P0RMAT(5X* SPECIMEN FAILED BEF0RE*E10 . 3*CYCLES*/) 

GO TO 1001 

1003 PRINT 39,FCYC 

39 PORMAT( 5X* THE SPECIFIED NUMBER OF LOAD CYCLES HAS^t 
1^ BEEN r®T*E11 .3) 

2C02 FORr«AT(/) 

1001 CONTINUE 
M=M/2 
N=N/2 

LOUNT=KOUNT-1 
PRINT 540 
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540 PORI^AT(/13X>^DSLTA-1^*■11]^TOTAL-I^ 18X-A"M6]?^DA/nr’< ISX^S^- 
1l6X*DB/rr">/) 

PRIZ'IT 5356, (FXIICY(I) ,PT!lC(l) ,PA 1 (I ) ,CPR( I ) ,FB 1 (I ) ,nPR(I 
1),I=1,L0UITT) 

536 PORI I AT( 6220.3 ) 

IP (TEST.IIS.O.) GO TO 1000 
EIID 

SUBROUTnra IpyER(S,M ,I>) 

c 

C THIS ROUTniE IRTERTS A 5X3 FATRIX 
C 

DIKaiSIOi: A( 3 , 3 ) ,AI ( 3 . 3) ,AG( 3, 3) , s( 3, 3 ) 
D=S(1,1)-^S(2.2)"S(3,5)+S(1,2)^'S(2,3^S(3,1)+S(1,3)* 
1S(2,1hS(3,2) 

2-fS( ip)!'S(3.l)^S(2,2)+S(l,2)^S(2,l)^‘S(3,3)+S(l,1> 
3S(3,2)’3(2,3)) 

AG(1,1}=S(2,2)*S(3,3)-S(2,3)^S(3,2) 
AG(1,2)=-S(2,1)^S(3,3)+S(2,3)^S(3,1) 
AG(1,3)=S(2,1)^:S(3,2)-S(2,2)^-S(3,1) 
AG(2,1)=-S(1,2)^S(3,3)+S(1,3)*S(3,2) 
AG(2,2)=3(l,lh'S(3,3)-S( 1,3)*S(3,1) 

AG(2,3)f-S( 1 ,1)>'S(3,2)+S(1,2)^S(3,1) 

AG(3,1)=S(1,2)' S(2,3)-S( 1,3hS(2,2) 

AG(3,2)=-S( 1,1)-S(2,3)+S(1,3)*S(2,1) 
AG(3,3)=3(1,1)>S(2,2)-S( 1,2>S(l,2) 

DO 1 1=1,3 
DO 1 J=1,3 
AI(I,J)=AG( J,I)/D 

1 CCUTIimE 
RETURN 
HID 

SUBICUTIrlE MULT(A,B,C) 

C 

C THIS ROUTINE I-ULTnilES TITO MATRICES 
C 

DIMENSION A(3,3),B( 3),C(3) 

DO 1 1=1 ,3 
SUM=0.0 
DC 2 r= 1 

2 SUi:=Glil+A(l,K)>B(K) 

c(i)=sun 

1 CONTIIIIE 
RETURN 
END 

SUBROUTEIE 3II Q( A,B,N,KS ) 

C SOLUTION OP THE LINEAR ALGEBRAIC STSTEII OP EQUATIONS 
C OP OHE PORI'I 

C AX=B 

C A - NXII MTRIX OP COEPPICIENTS 

C B - VECTOR FOR THE RIGKTHAND-SIDE OF THE SYSTEM OF 

C EQUATIONS 
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N- NUMBER 0? EQUA.TIONS AND VARIABLES 
KS - OUTPUT DIGIT 

0 POR A NORI'iAL SOLUTION 

1 FOR A SINGULAR SET OF EQUATIONS 
DIMENSION A(1),B(1) 

T0L=0.0 
KS=0 
JJ= -N 

DO 65 J=1,N 
JY=J+ 1 
JJ=JJ+N + 1 
BIGA=0.0 
IT=JJ-J 
DO 50 I=J,N 
IJ=IT+I 

IF(ABS(BIGA)-ABS(A(IJ))) 20,50,50 
20 EIGA=A(U) 

IMAX= I 
30 CONTINUE 

IF(ABS(BIGA)-TOL) 35,55 ,40 
35 IS=1 
RETU.^T 

40 I1=J+N-(J-2) 

IT=IMAX-J 
DO 50 K=:J,N 
11=1 1+N 
I2=Il4lT 
SAVB=A(I1) 

A(I1)=A(I2) 

A(I2)=SA’C 
50 A(ll)=A(Il)/ErGA 
SAVE=B(IT>'IAX) 

B(irAX)=B(J) 

B( J)=SAVE/BIGA 
IF (J-N) 55,70,55 
55 IQS=N>^(J-1) 

DO 65 K=JY,N 

IXJ=IQS+IX 

IT=J-IX 

DO 60 JX=JY,N 

IXJX=N^(JX-1 )+IX 

JJX=IXJX+IT 

60 A( I X JX ) =A ( IX JX ) - (A ( i:{J ) ■> A( J JX ) ) 

65 B(IX)=B(IX)-(B(J) A(IXJ)) 

70 NY=N-1 
IT=IJ^N 

DO 80 J=1,NY 

IA=IT-J 

IB=N-J 

ic=i: 

DO 80 N=1,J 
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B(IB)=B(IB)-A(IA)^B(IC) 

IA=IA-N 
80 IC=IC-1 
RETURN 
END 

SUBROUTINE PORM(N,M) 

THIS ROUTEIE FORI-IS THE COEPHCEINT MATRICES USED TO SOLVE 
FOR ELASTIC BITERIAT^INAR STRESSES 

ETMEJISION THETA ( 10) 

C0MM3N/XTjIMIT/LC( 2,51 ),DB(6, 51 ) ,EP 

COr-a^50N/ROT/Z9(lOOOO),D( 100 ) ,mLF 

COM'DN/TOr /E3, T2 , V 1 , 3-IX , 3'IY , G , CONS , Q , A1 , SCON 

COia-DN/BOT /E4, T4 , V4 , SCX , SCY , GB , OONB , Q 1 , GC , Ep 

COI'MOII/ADHES/TAU ,GAD 

CORCN/BOND/P , PI , P2, XKF ,XKUNS ,XKSTIF 

C0MM0N/B0ND2/NE ,NL ,NT ,XD( 100 ) ,TC( 100 ) ,XA( 100 ) , YA( 100 ) , 
INOP(IOO) 

COMTON/CTOL/TOL.NC ,m ,TX ,TY ,NBC( 100),LBC 
COMPLEX Z,CI,F1,G2 
EXTERNAL F1,G2 
CI=CMPLX(0. ,1.0) 

A=A1 

0=M/2 

P=N/2 

L=P 

MX=NC*10l 

K=0 

ASSEMBLE RIGHTHAND SIDE VECTOR 
PRINT 57 

57 FORMATC /- node NUI'IBER*8X*X-OOR^8X->'-Y-COR*/) 

PRIIff 56, (I ,DC(1,I) ,DC(2,I) ,1= 1,MX) 

56 FORMAT(I10,2E15.5) 

DO 101 1=1, K 
X=DC(1,I) 

Y=DC(2,I) 

Z=X+CI*Y 

CALL REMIU(DU,DV,Z) 

CALL REIOT(BU,BV,Z) 

D(I)=IU-BU 
D(K+I) = DV-F\T 
DO 102 J=1,MX 

CALLXINTG(Z,J,A,C1,C2,C3,C4,F1,1 ) 

CALL nNTG(Z,J,A,B1,B2,B3,B4,P1,2 ) 

802 CONTI.TOE 
11=1 
J1=J 

IV1=(J1-1 + 
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Z9(IV1)=C1+B1 

IF( I.BQ.J)Z9(IV1)=Z9(IV1)+TAD/GAD 
12=1 
J2=J+L 

IV2=(J2-1 ^N+I2 

Z9(IV2)=C2+B2 

I3=I+K 

J3=J 

IV3=(J5-1 H"+I3 
Z9(IV3)=C3+B3 
14=1 + K 
J4=J+L 

IV4=(J4-1 )*N+I4 
Z9(IV4)=C4+B4 

IF(I4.BQ .J4)Z9(IV4)=Z9(IV4)+TAD/GAD 
102 CONTINUE 
101 CONTINUE 
104 CONTINUE 
20 CONTINUE 
RETURN 
END 

SUBROUTINE REMTU(U,V,Z) 

CALCULATE DISPLACEIIENTS IN THE TOP SHEET DUE TO A REMOTE 
STRESS 

OOMMCN /TO P/E3, T2 , VI , 3DC , aiY , G , OONS , Q , A1 , SCON 
C0I4PLEX CI,Z,ZI,ZB,CXSR,D,IiPHI ,PHI ,DOMEGA 
CI=CMPLX(0. ,1.) 

ZB=ODNJG(Z) 

D=CXSR(Z,A1) 

DPHI= . 5 *SM Y*D- . 25*Z* ( SHY-3DC ) 

PH I = . 5 *SM Y* Z/D- . 25 *( SMY-^IX ) 
D0MEGA=.5*Sm’^CXSR(ZB,Al)+.25 *(SMY-3U0*ZB 
D= ( Q ><^DPH I -DOME GA- ( Z - ZB ) * C ON JG ( PHI ) ) / ( 2 . -»G ) 

U=REAL(D) 

V=AIT4AG(D) 

RETURN 

END 

SUBROUTINE REI4BU(U,V, Z) 

CALCUIATES DISPLACH'IENTS IN THE BOTTOM SHEET DUE TO A 
RII40TE STRESS 

COMMCfl /BO T /E4, T4 , V4 , SOX , SCY , GB , CONB , Q 1 , GC , E5 

COMPIEX Z 

X=RIAL(Z) 

Y=AIMAG(Z) 

S12=-V4/E5 

U= ( SCX^5+S 12 *SCY }*X 
V= ( S12*SCX+30Y/E4 >Y 
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RETURN 

END 

COMRLEX FUNCTION F1(Z,Z0) 
EVALUATE INTEGRAL OF B(Z,Z0) 


COMILEX Z , ZO , ZB , ZOB ,XB ,XD ,CXSR ,CI , D,B,R 

COMPLEX 

DIM2NSI0N N(30O) 

C0IWDN/T0P/E3, T2,V1,35X,aiY,G, CONS , Q , A1 , SCON 
COMMDN/XK/KDOUNT 

COUNTER LIMITS INTEGRATION INTERVAL TO 50 STEPS 
USE COMMCN TO ZERO COUEITER 

XA(Z,Z0)=CL0G((Z^Z0-A**2-GXSR(Z,A)*CXSR(Z0,A)y(Z*Z0+A 
2 ** 2 - 

2CXSR( Z ,A)*CXSR( Z0,A) )* (-1 . )) 

XB( Z,Z0 ) =CL0G( (Z^Z0-A^^*2-CXSR( Z,A)*CXSR( Z0,A) )*( Z+ZO )* 
1 * 2 / ( 

2 (Z* ZOfA* 2- GXS R ( Z , A) CXSR( ZO , A) ) ^ ( ZO-Z) ** 2 ) ) 

XC ( Z , ZO ) = ( Z-ZO* CXSR( Z , A) /CXSR ( ZO , A) ) * ( ZOB-ZO ) / ( Z* ^ 2- 
1ZO**2) 

G1=G 

ZB=OONJG(Z) 

Z0B=C0NJG(Z0) 

CI = CMPLX(0. ,1. ) 

SAVE PREVIOUS VALUES OF XB ML CHECK PATH FOR BRANCH CUT 
IF BRANCH IS CROSSED ADD PROPER IMAGINARY TERM TO KEEP 
DISFLACEI'IENTS CONTINUOUS 

KC0UI^T=:KC0UNT+1 

K=KOOUNT 

A^AI 

BI = .5*(XB(Z,Z0)-Q»^XB(Z,Z0B) )+XC(Z,Z0) 
X=REAL(XB(Z,ZOB)) 

Y=AIIAG(XB(Z,ZOB)) 

1=1 

IF(X.LT.0.AND.Y.GE.0)I=2 
IF(X.LT .0 .AND.Y.Iff .0)1=3 
D=XB(Z,ZOB) 

IF(K.LE.2) GO TO 100 

IF(N(K-2).B3 .2.AND.I.IQ.3)D=D-6.283*CI 
IF(N(K-2).]33.3.AND.I.B3.2)D=D+6.283*CI 
100 CONTDIUE 
N(K)=I 
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H1 = -Q'XC(Z,ZO) 

H=-Q'<REA1(XA(Z,Z0)) 

H2= . 5 *( Q ^ C OI'UG ( I) ) ) 

R=H+H1+H2 

F 1 =5l+XC( ZB , ZO ) + 2*( Z^: ZO-ZOBi^" ZB ) /(ZB^- ^ 2- 30^-<^2 )+ ( Z -ZB ) « E ( 
1ZB,ZO) 

G1 = 5. 75^10*^ 6 
F1=F1^-C0IIS 
D1«=REAI(D) 

D2=Airi/iG(D) 

RETURN 
EITD 

JUBROUTIIIE XINTG(Z,I,A,S11 ,S12, S21 , S22 ,F1, IC) 

mis ROUTINE INTEGRATES ffiE SUI-! AND DIFFERENCE OP !R70 
COr-E LEX FUNCTIONS, F1 AND F2, USED 'JD EVALUATE THE GREENS 
FUNCTIONS FOR DEPLACEXEI^TS. 

COW'DN/CTOL/TOL , NC , M , TX , TV , NBC ( 100 ) , LBC 
CO]CnON/IX/KCOUNT 

COIffiTiEX F1,F2,Z,Z0,CI ,A1,B1,Z0B 
COI-iroiT/NLIKIT/DC ( 2 , 51 ) , D3 ( 6 , 51 ) , EF 
EXTIRNAL FI 
CI=a-I?LX(0. ,1.) 

X=REAI.( Z) 

311=0.0 
312=0.0 
321=0.0 
322=0.0 
Y=AU*AG(Z) 

C 

C IF u:ds=i tii-fi z z-ies outside iiwegration patch 

C AITDNI 10 1 NTS ARE USED HI A SBIGLE INTEGRATION. 

c IF d:dex=2 bien msN z lies ’nitiiin mE integration patch 

C A17D TNO INTEGRATIONS ARE HADE WITH K2 POINTS IN EACH 
C LNTSdlATIOI: . 

C ICOUNT CONTROLS THE LOOPS ON THE INTEGRATION. IF ICOUNT 
C EQUALS 0 CONTBIUE GTHExIWISE DO SECOND HITSGRA'EION . 

C 

:<ST/aIT=DB(l,l) 

:ai’INAL=D3(2,l) 

IF ( X .IE .X5 TART . OR .X . GE .XFINAL) INDE:':= 1 
I? ( X . GT . XS TART . AND . X . IT . XF I NAL) INDEX= 2 
IC0UIIT=0 
Id = 5 
K1 = 5 

IFONDEX.IQ.1 )2,5 
2 .&=}START 
.■if=:vPii:aj. 

GO TO 4 
5 yi=XSTc7RT 
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XU=X-T0L 

icouNT=icourmi 
GO TO 4 

5 XU=XFHIAL 
.a=X+TOL 
IGOUITT=ICOUITT+1 

4 COIITIIU3 

IP(:CPII';AL.B3 .A)XFUIAi=A-TOL 

iK=(:{u-xL)/(3.*i:) 

M=K+1 

K) 7 L= 1 ,n 

XX=(XU-XL)^(L-1 ) 

XD = XX/K+XL 

X1=AT/2. 

X2=L/2 

.0R.L.B3.M)16, IS 

16 XF=1 

GO TO 17 

18 IP(iffiS(X1-X2).IT..0001 )XF=4. 

IP(AES(X1-X2).GT..0001 )XP=2. 

17 COUTHIUE 

zo=xofcr 0.0 

y0=AII'IAG(Z0) 

ZOB=COnJG(ZO) 

KC0mT=0.0 

CALL YI1ITG(Z,Z0,I,A,T11 ,T12, T21 , T22 ,F1, IC) 

201 PORKATCA Y-STRIF at XD=*E10.5^T1,T2,T3,T4=^^4E11 .5) 
C11=XF^DX^T11 
C12=XF*DX-T12 
C21=CCF^;-DX^T21 
C22=XF^D]^T22 
311=3114011 
312=312+012 
321 =321 4021 
322=32240 22 

7 con T HUE 

IF(IC0U1IT-1 )6, 5 ,6 

6 coNTnros 

RETUPJI 

EtID 

3UBR0UTIIS YI1ITG(Z,Z0,I,A,S11 ,312, 321 , 322 ,F1 , IC) 

THIS ROUTIIIE IIITEGRATE3 THE 3UM AHD DIFFERENCE OF TWO 
COMPLEX FUITCTI01I3, FI AND F2, USED TO EVALUATE THE GREENS 
FUNCTIONS FOR DISPLACEI'NSNTS. 


COMHjEX H,A2,B2 

COIFON/CrARI'I/S1,S2,C11 ,C 12 , C22 ,C21 , PI , P2, Q4, Q2 
COMPLEX Z1 ,Z2,W1 ,\/2,G,GB, 31,32, C11 ,C12, C22 ,C21, PI ,P2, Q 
14,Q2 
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OOMMON/XLIMIT/IC ( 2 , 51 ) , IB( 6 , 51 ) , EF 
COMMCJT/CTOL/TOL ,KC ,NR ,TX .TT ,KBC( lOO ,IBC 
COMPLEX E1 ,P2,Z,Z0 ,CI ,A1,B1,ZOB 
COI-IMai/RDT/Z9( 10000 ),D( 100 ), HALF 
CI=CF-PLX(0. ,1 .) 

XD=RSAL(3D) 

X=REAL(Z) 

311=0.0 

312=0.0 

321=0.0 

322=0.0 

Y=A]I'IAG(Z) 

B1=DB(3,I) 

B2=DB(4,I) 

C1=0B(5,I) 

C2=DB(6,I) 

Y1=YYX(X0,B1,B2) 

Y2=YYX(X0,C1,C2) 

IC0UNT=0 

K1=3 

K1 = 5 

XL=Y1 

XU=Y2 

11 = 1 

IF(n.B3.l)G0 TO 4 
C 

C lOGIC TO STATEIffiKT POUR IS FOR THE NARROW STRIP THAT 
C OONTAINS OHE SHTGULARITY 

C THESE STATEMENTS CAN BE USED TO INTEGRATE THIS STRIP WITH 
C ADDITIONAL LOGIC. CONSIDER X CONSTANT OVER THIS NARROW 
C STRIP 
C 

IF ( Y. IE . Y 1 . OR . Y . (2) . Y2 )INDEX= 1 
IF(Y.GT .Y1 .AND .Y.IT .Y2)INDEX=2 
IF(irnEX.B3.l)2,3 

2 :CL=Y1 
XU=Y2 
GO TD 4 

3 :CL=Y1 
XTJ=Y-TOL 
IC0UNT=IC0UNT+1 
GO TO 4 

5 XL=Y+TOL 
:0J=Y2 
K1=5 

IC0UNT=IC0U1IT+1 

4 CONTINUE 
K=K1 

DX=(XU-XL)/(3.*K) 

M=K+1 

DO 1 L= 1 ,M 
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XX=(XU-XL)*(lr-1 ) 

YO=XX/K+XL 

AT=L 

X1=AT/2. 

X2=L/2 

IF(L.B3.1 .0R.L.H3.M)16,18 
16 XF=1 

GO TO 17 

18 IP(ABS(X1-X2).IIP..0001)XF=4. 

IF( ABS (X1-X2) .GT ..0001 )XF=2. 

17 CONTETUE 

ZO=XQfCI^O 
ZOB=CaiJG(ZO) 

IC=1 CRACfZED lETAL SHEET 
10=2 ORTHOTROHC SOHD SHEET 

IF(IC.KE.2)30,31 

30 OONTIIIUE 
A1=F1(Z,Z0) 

B1=F1(Z,Z0B) 

T11=REAL(A1+B1) 

T12=REAL(CI*(A1-B1)) 
T21=AIMAG(A1+B1) 
T22=AiraG(CI*(A1-Bl)) 

GO TO 33 

31 OONTIIKJE 
Z1=X+S1*Y 
Z2=X+S2*Y 
W1=X0 + S1*Y0 
W2=X0 + S2*-Y0 
Ai=G( zi ,vn) 

B1=<J( Z2,V/2) 

A2=H(Z1 ,W1 ) 

B2=H(Z2,U2) 

T 1 1 =2 . *REAL(F 1*0 11 *A 1+P 2*0 21 *B 1 ) 
T12 = 2.*REAL(P1*C 12*A2+P2*C22*B2) 

T21 =2.*REAL(Q4*C 11*A HQ2*C21*B1 ) 

T22 = 2 . *REAL (Q4^C 1 2 *A 2-fO 2*0 22*B 2 ) 
33 OOIITETUE 

IF(I0.HE.3)60,61 

60 FS=1. 

FR=1. 

GO TO 62 

61 OOIITINUE 
FS=D(I) 

FR=D(I+1TALF) 

62 OOllTIIJUE 
P11=XF*FS*DX^T11 
P12=XF^<FRXDX*T12 
P21=XF*FS^DX*T21 
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?22=XF*FR^DX^T22 
S11=S11+?11 
312=312+1-12 
321=3 21 +?21 
322=3 22+? 22 
1 OOIITIilUE 

IF(IC0UMT-1 )6,5,6 
6 CONTIIIUE 
RETUH'I 
HID 

COMiliEX FUITCTION CXSR(Z,B) 

C 

C THIS FUNCTION TAKES THE SQUARE ROOT OF Z*^2-B**2 

C IT ONLY RETURNS OHE FIRST ROOT, THE SECOND ROOT CMl BE 

C FOUND BY ADDING H /2 TO THE TRIP ARGUI'IENT 
C 

COMPLEX Z,T1,T2,CI 

ci=ci>:pix(o. ,1.) 

T1=Z-B 

T2=Z+B 

X1=REAL(T1) 

Y1=AIM/JG(T1) 

:C2=REAL(T2) 

Y2=AII'iAG(T2) 

A1=ATAII2(Y1,X1 } 

A2=ATAN2(Y2,X2) 

S1=(REAL( Tl)^f ^2+AIMAG<T1>*2>* .5 
S2=( REAL ( T2h " 2+A IMAG( T2 )^ * 2 .5 

3R=(S1*S2)^' -> .5 
AIIG=A1+A2 

CXSR=SR>' (C0S(ANG/2 .0)+CI^^SIN( ANG/2 . )) 

RETURN 

END 

COMPLEX FUl'ICTION XK(Z,ZO) 

C 

C GREEN'S FUNCTION FOR STRESS DITENSITY 
C 

CCMMCZ: / ton/e:, T2, VI , SI'IX , 3-lY , G, CONS , Q , A1 , SCON 

CON'l'IE}: Z, ZO , ZB , ZOB , CXSR , D, F 

A=A1 

ZB=ODNJG(Z) 

ZOB=CCNJG(ZO) 

D=(ZO' Z0B-2.-xZ0'->:2+A^^2)/((Z0*’>'2-A^^ 2)»CXSR( ZO,A) ) 

F=Q ^CI{SR(ZOB,A)/(ZOBx--^2-A^*2) 

:G:=2^-A^" .::-»'-(F+D)/(6.283»*( 1.4Q)-s<T2) 

RETURN 

END 

FUNCTION YY(T) 

C 

C T IS A DUI'KY PARAMETER FOR XO 
C CALCUIATE YO FOR A GF/Ei: XO 
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COIMOIT/TOI /E3, T2, V1 , 3IX , 3-1? , G, CX3NS , Q , A1 , SCON 
COinvlOiT/BOND/F , P1 , P2, HCP ,XKU1IS .XI^STIFF 
I?(T/A1.LE.O)2, 1 
1 CONTIMUE 

YY=I» ( ( 1- (T /A 1 n )■« ^ ( 1 /P2 )) ^-A 1 

GO TO 3 
2 yy=5"'<A1 

3 CONTI HUE 
RETUlUJ 
END 

cor-iPiEx fui7Ctio:j b ( z , zo ) 

COI MON /TO P /E3 , T 2 , V 1 , 3'IX , 3IY , G , CONS , Q , A1 , SCON 
COI :i TEX X 1 , X2 , X3 , X4 , X5 , CXS R , ZOB , 20 , Z , X6 , X7, X8, X9 
01=0 
A=A1 

ZOB = COIIJG(ZO) 

X1=CXSR( Z,A) 

X2=CXSR(Z0,A) 

X5=1./(Z-''- ''2-20>'^2) 

X4=1./(Z-^! 2-Z0B->! -2) 

X5=1./(Z-Z0)^ --'2 
X6=1./(Z+Z0)>^-"2 

X7=X 1 *( ( -4 *Z0 * X3 )+ 2 . ’'O 1 *ZOB* X4+( Z - ZOB ) ^ X 5- (Z+ ZO B ) X6 ) 
X8=( ZO-ZOB ) > ( ( A^ > 2- Z^ZO X5- (A* 2+Z^ ZO )*X6) /X2+2. *X 2^ 
1X3*Z 

X9= -2 . ^ 1 *Z^- C XSR ( Z OB , A) ^ X4 
B= CC7+X8+X 9 ) / (2 . 0*X1 ) 

RETURN 

EITD 

SUBROUTINE CPAR 

CALCUIATE lARAHTERS FOR COM l-' LEX VARIABLES IN ORTHOTROPIC 
ANALYSIS. 


COMilEX Cl ,31,S1B,S2,S2B,C11 ,C12, G21 , C22 ,I'1,P2, Q4, Q2,D 
1,U1,U2 
COIHLEX CC 

C0I3»10N/B0T/E4, T4, V4, SCX , SCY , GB , CONB , Q1 , GC , E5 
COI'uION/C: ARi:/31 , S2,C 11 ,C12, C22 ,C21 , 1 1 , P2, Q4, Q2 
CI = CIRLX(0.0, 1.0) 

EX=E5 

EY=E4 

FX=V4 

GXY=GC 

B=X.:/GZY-2*TA 

C=EX/EY 

CC=CM1’LX(C,0.0) 

VY=vx) EY/s;: 

D=CS0RT(B-'''2-4.*GC) 

U1=-.5»^(B-D) 
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U2=-.5*(B+D) 

S1=CSQRT(U1) 

IF(AII^AG(S1) .Iff .0 )S 1=-S1 
S2=CSQRT(U2) 

IF ( Al'IAG (S2 ) . Iff .0 )S 2=-S2 
S1B=C0NJG(S1) 

S2B=C0NJG(S2) 

IF(F^ 2-4. -K3) 1,2,3 

1 CONTINUE 
PRINT 4 

4 FORKAT(* ROOTS ARE COMPLEX U1^+IB, U2=-AfIB*/) 

GO T05 

2 CON T HUE 
PR HIT 6 

6 FO?J.U\T(^ ROOTS ARE PURE IMAGINARY AND EQUAL*/) 

GO T05 

3 B1=AIMAG(S1) 

D1=;AIMAG(S2) 

P1=(S1**2-VX)/EX 
P2-(S2**2-YX)/EX 
Q2=( -VY* S2+1./S2)/Ey 
Q4=(-VY*S1+1./S1)/ET 

100 F0RHAT(4E12.3) 

C0MP1=1 ./(12.566*T4»(B1**2-D1**2)) 

C 1 1 =( D 1 ** 2*VY+ 1 )* B 1*C0MP 1 
C12=(VX+D1**2)*C0MP1 
C21 =-D1*( VY*B1** 2+ 1 )*C0MP1 
C 22 = - (B 1 ** 2+VX ) *C OMP 1 
GO TO 7 

5 HUNT 8 

8 FORI'IAT(/* ERROR-ERROR QRTHOTROPIC ANALYSIS IS NOT* 

1* DEFINED ■></) 

7 CONTINUE 
RETUHT 
END 

SUBROUTINE GRID 

GENERATE MESH USED TO DISCRETIZE INTERLAMINAR STRESSES. 

COMMON /B0ND2A"E,NL ,NT ,XC( 100 ) , YC(100 ) ,XA( 100 ) , YA( 100 ) , 
INOP (100) 

OOMMCN /CTO L/TOL , NC , NR , TX , TY , NBC ( 100 ), LBC 
COraiON/TOP/E3, T2, VI , »IX , S^IY ,G, OONS ,Q , A1 , SCON 
COMMa: /:CLIM IT/DC ( 2 , 51 ) , IB ( 6 , 51 ) , EP 
COMMON /BOND/F , PI , P2, XKF ,XKUNS ,XKSTIF 
P=P1 

TXX=2.*T0L 

DX=TX/1TC 

DY=TY/NR 

LBC=0 

MX=NCXNR 
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DO 11 1=1, M 
11 NBC(I)=0 

DO 100 J=1,NR 
DO 100 1=1, NO 
MT=NC*(J-1 )+I 
AG=A1+DY*(J-1 )+DY/2. 

YI = F*A1 

BC=YI+DY^(J-1 )+DY/2. 

DO ( 1 ,MT )=DX/2 .+( 1-1 )^DX 
jt=DC(l,MT) 

Y=YYD(X,BC ,AC ,P) 

DC(2,MT)=Y 
DB(1,MT)=DX''(I-1 ) 

DB(2,KT)=DX^‘I 
DB(3,I®)=DY*(J-1 )+YI 
DB(4,MT)=D^(J-1 )+A1 
DB(5,M)=IY*J+YI 
DB(6,ME)=Dr*J+A1 
100 COITTIIIUE 
110 P0RMAT(I10,2E15.3) 

120 F0RMAT(I10,6E12 .3) 

RETURN 

END 

SUBROUTINE XING (SUM, XNCY, IAS , DPS ) 

THIS RDUTHIE INCREMENTS THE CRACK LENGTH Al'ID DEBOND SHAPE 
PCR NCY CYLES 

CCM-DN/XLIMIT /DC ( 2 , 51 ) , DB( 6 , 51 ) , EF 
CO MMO N /CT 0 L/TO L , NC , NR , TX , TY , NBC ( 100 ) , LBC 
COMMON/BOND/F , ?1 , P2, XKF ,XKUNS .XKSTIFF 
COMMON /TO r /E3 , T 2 , V 1 , aiX , 31Y , G , CONS , Q , A 1 , SCO N 
COMM® /DO T /S4 , T4 , V4 , sex , 9CY , GB , CCNB , Q 1 , GC , E5 
COMffiN/XY/XD ,YD 

DIMENSION YN( 20 ) ,XII( 20 ) ,ET ( 20 , 2 ) ,CI ( 20 ) ,D( 20 ) ,CR( 2, 2) 

B1=A1*F 

R=.01 

DA=3. 22E-14 ^XKSTIFF**3 .38 
CFAC=1.79E-14 
CFAC=3. 36E-14 
DA=CFAC^XKSTIFF^<*3. 38 
m=DA/(( 1. -R)*56000 .-XKSTIFF) 

DF=3. 158E-05*SUM-x--''3.616 

DAS=IA 

DFS=DF 

DETERMDIE HOW MANY CYCLES REQUIRED FOR EITHER A CRACK OR 
DEBOND EXTENSION OF .1 INCHES. THEN USE SMALLEST VALUE 
AS THE INCREMENT OF APPLIED LOAD CYCLES 


XNCRACK=.10/DA 
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XNB01'3D=.10/IIF 

XNB0ND=.20/IF 

ja;CY=XNCRACK 

IF(XNGRACK.G!P .XNBOirD)XNCY=XNBOIJD 
DA=X1ICY»^DA 
A2=^1+DA 
DF=DF*XNCY 
B2=B1+DF 
?1 = 2 
A1=A2 
F=B2/A2 
RBTURIJ 
EIJD 

SUBBDUTIME VERI ( Z , SXX , SXY , SXXY , SYX , SYT , SJfXY ,MX , KTYPE ) 
C 

C INTEGRATE GREENS FUNCTIONS FOR STRESSES 
C 

COMMOIT/CTO L/TO L ,NC , NR , TX , TY , NBC ( 100 ) , IBC 
COMPLEX Z 

z= z 
MX=r-K 
J=J 
SXX=0. 

SXY=0. 

SXXY=0. 

SYX=0. 

3YXY=0. 

SYY=0.0 
DO 1 J= 1 ,I4X 

CALL VXB'ITG ( J , Z , SI , S2, S3, S4, S5, S6, MX , KTYPE ) 

C 

C FIRST INDEX INDICATES LOAD DIRECTION, SECOND STRESS 
C DIRECTION 
C 

3XX=S1+SXX 
SXY=S3+SXY 
SXXY=S5+SXXY 
SYX=S2+SYX 
SyY=34+SYY 
SYXY=S6+SYXY 
1 CCNTB'US 
RETURN 
END 

SUBROUTDIE VXINTG(l , Z , S1 , S2, S3, S4, 35, S6, 14X , KTYPE) 

C 

C DETERI'IINE STRESSES IN ADHHIENDS DUE TO INTERLAMINAR 
C STRESSES. 

C 

COMON/XLIMIT /DC ( 2 , 51) , IB ( 6 , 51) , EP 

C0FJ-r)N/R0T/Z9( 10000 ) ,D( 100 ) ,NALF 

COKMaj /TO P/E3 , T2 , V 1 , 3DC , 3'IY , G , CONS , Q , A 1 , SCON 
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C 

C 

c 

G 

c 

c 

c 


COMILj]}: Z ,Z3,20 

CO! iMa: /CTC iVto l , iic , m , ®c , w , nbc ( ioO , isc 
COI'.Il'L.'vX Cl 
GOIilLBX i' 
ii=z 

':=RZAL(Z) 

31=0 


32=0 

35=0 

34=0 

35=0 

36=0 

ZB=CO:iJG( z ) 


Go f G7 f G3 


IF IIIDIX=1 TIL3N 2 TZE OUTSIDE THE lUTEGRATIOlI BATCH 
AIH) i:i lOIKTG AROUSED IN A SINGLE INTEGRATION 
IF IIDEX=2 THEN Z I lES NTTHIi: THE INTEGRATION lATCH AND 
TNC INTEGRATIOi:S ARE IIEADE EACH WITH K2 INTEGRATION POINT 
IF IC0UNT=0 TIiSi: OONTINUE OTHERWISS SECOND IKTEGllATION 


.'3TART = DB(1 ,1) 

XFI1]AL=I!D(2,I) 

IF ( Z . IE . :(S TART . OR . X . G3 . ZFI NAL) INDEX= 1 
IF(Z.GT ECSTART . AND. X. Iff .XFINA].)IND1C{=2 
K1=5 
K1=5 

icou:iT=o 

IF(INDEX.BQ .1 )2, 3 

2 XL=X3TART 
yU=XFINA: 

GO TO 4 

3 Xr-=XSTART 
.aj=x-TC'i 
IC0UNT=IC0UNT+1 
GO TO 4 

5 :<U=XFINiil. 

/Jj=^k+T0_j 

IC0U]3=IGC'U]ITl-1 

4 CONTINUE 

IF (IT INAL. EQ .A1 )ZFINAL=A1- TOL 
N=K1 

dx=(xu-:c)/(3.'h. ) 

K=X+1 

ci=a:iTX(o. ,1 .0) 

DC 7 L=1,H 

:C':=(xu-:t)^'(l-i ) 

.:o=.a/K+z-. 

•TO=yY(xo) 

ZC = XC+CI''YO 
/,03=coi:jg(zo) 

AT=T 
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X1=AT/2. 

X2=L/2 

IF(L.iD3 .1 .OR .L.S3 .M) 16, 18 

16 ZP=1 

GO TO 17 

18 IP ( i'iB S (X 1- X2 ) . IT . . 00 01 )X?= 4 
IP(ABS(X1- X2) .GT..0001 )XP=2 

17 CONTETjE 

CALL VYINTG( Z ,Z0 , 1 , TXX , TXY , TZX , TYT , TXYX , TXYY ,mPi: ) 

S1=XP^E{> TXX-f SI 

S2=XF->'IEv-' TXY+32 

33=XF->-IK^TmS3 

34=XF->:-rK>'TYY+34 

S5=XF>^IK'- T:YX+3'-. 

36=XF-IK*TXYY+S6 
7 CONTIinjE 

IF(IG0U11T-1 )6,5 ,6 
6 COUTINUE 
RET URL 
SID 

SUBHOUTDTE VYIIITG ( Z,ZO , 1 , 37, 38, 35, S4, S3, S5,KTYFE) 
INTEGRATE GREENS FUNCTIONS FOR STRESSES 


COHMOll/Ch'ARI'l/S 1 , 32, C11 ,C 12 , C22 . C21 , PI , F2, Q4, Q2 
COHmi /XL IHIT /DC ( 2 , 51 ) , DB ( 6 , 31 ) , IF 

CONHEX Z1,Z2,W1,W2, GB,S1,S2,C11 ,C 12, C22 ,C21 , II ,P2, Q 


14, Q2 

COXI'Di:/TC P /E5 , T2 , V 1 . 3-IX , S^Y , G , CONS , Q , A 1 , SCO II 
COMI'Di:/ROT/Z9(lOOOO),D( 100 ),MLF 

COMlXEX Z,Z3,ZO,ZOB,B,IBZ,Bl,r!B1B,B1B,DB1,G5,G5,G7, G8 
COMPLEX CI,G1A 

COIIPXEX G 1,G2,G5,G4,H1,H2,W1B,W2B 
CCMI-iai /CTO L/TO L , NC , NR , TX , OY , NBC ( 100 ) , IBC 
CCXIXEX P 


X=REAI,( Z) 
Y=AIMAG(Z) 

:=z 


Q1=Q 

ZB=CONJG(Z) 

37=0. 


33=0.0 

;3=o 

34=0 

33=0 


36=0 

I'E{=ra*NC 


DB GIVES PARAIETER3 


FOR UPPER AND IDV/ER GRID 30LUroARIE3. 


B1=DB(3,I) 
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B2=DB(4,I) 

C1=DB(5,I) 

C2=DB(6,l) 

XO=REAL(ZO) 

Y2=YYX(XO,C1,C2) 

Y1=YYX(:-3D,B1,B2) 

E1 = 5 
El = 3 
JCL=Y1 
X(J=Y2 
E=E1 
C 

C SEE YIIITG FOR LOGIC FOR INTEGRATION OF NARROV; STRIP THAT 
C OONTAH'IS THE SINGULARITY. THAT STRIP IS NEGLECTED HERE 
C 

DX=(XU-XL)/(3.-^) 

M=E+1 

CI = CK?IX(0. ,1.0 ) 

DO 7 L=1,M 
XX=(XU-XL)^-(L-1 ) 

YO=XX/K+XL 

ZO=XO+CPYO 

ZOB=OONJG(ZO) 

AT=L 

X1=AT/2. 

X2=L/2 

IF(L.E3 .1 .OR.L.KJ .M) 16, 18 

16 :<F=1 

GO TO 17 

18 IF(ABS(X1-X2).IT..0001 )XF=4 
IF(ABS(X1-X2).GT..0001 )XF=2 

17 CONTINUE 
IFCKTXPE.iQ .1 ) 30, 31 

30 OONTIIRJS 
B1=B(Z,Z0) 

B1B=B( Z,ZOB) 

DB1=DBZ(Z,Z0) 

DB1B=DBZ( Z.ZOB) 

G1 = -4nZ0/(Z^*2-23**2)+Z0B/(Z^*2-Z0B**2)) 
G1=REAL(G1) 

G1A^-4*(Z0/(Z^->=2-Z0* >«-2)-ZDB/(Z->'*2- Z0B**2)) 
G1A;=REAL(CI^G1A) 

G2=2^ 1*ZOB/(Z^^'2-ZOB**2) 

G2=G2- ((Z0BfZB)/(Z+Z0)^^*-2+(Z0B-ZB)/(Z-Z0)**2) 
G5=2-!^ 1*Z0/(Z--^2- Z0**2) 

G3=G 3- ( (ZO+ZB ) /(Z+ ZOB ) 2+( ZO-ZB ) / (Z-ZOB ) 2 ) 

G5=B1+B1B 
G6=CI^* (B1-B1B) 

G7=(ZB-Z)-x (DB1+DB1B) 

G8=CI- (ZB-Z)’' (DB1-DB1B) 

T1S=REAL(G1- G2-G3) 
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T2S=REAL(G1A-CI^ (G2-G3)) 
T3S=REAL(G14G24G3) 

T4 S=RE AL( G 1A+ C ( G2- G3 ) ) 
T5S=AII-IAG(Gl4G2-fG3) 

T63=AIMAG( G1 AfC I-*( G2- G3 )) 

T11 = (-REAL(2^5-G7HT1S)^^SC0N 
T22 = ( -REAL ( 2*G 6- G8)+T 23) * S CON 
T3 5= ( -RE AL ( 2^K} 5-tO 7 )+T 3 S) ^ S CON 
T44=( -REAL ( 2*G 64G8 )+T 4 S) SCON 
T55=(AII'IAG(-G7)+T5S)*SCON 
T66 = ( Al-IAG ( -G8 )+ T 6 S)* SCON 

CHAIJGE SIGNS TO ACCOUNT FORNEGATIVE BODY FORCES 


T11=-T11 
T22=-T22 
T33=-T33 
T44=-T44 
T55=-T55 
T66=-T66 
GO TO 32 

31 CONTINUE 
Z1=X+S1*Y 
Z2=X+S2*Y 
V/1=XO+S1^YO 
V/2=XO+S2'YO 
V/1B=CONJG(vn) 

1V2B=C0IIJG(V^2) 

G1 = 2*( W 1 / (Z 1 2- v;i** 2 )+Vn B/ (Z 1 ** 2- * 2 )) 

G2=2*(U2/(Z2^ 2- V/2**2)+W2B/ (Z2^ 2-W2B^^ 2)) 
H1 = 2*(V/1/(Z1-'-'' 2- V;i** 2)-V/1B/(Z1*^ 2-W1S^*2)) 
112= 2*( 1, 2/ (Z 2^ 2- V2** 2 )-W2B/ (Z 2** 2- W2B>‘ 2 )) 
T11=2-’REAL(S1**2^11*G1+S2*^ 2^:0 21 ^2) 

T22 = 2*REA L ( Cl * ( S1 ** 2*C 1 2 *H 1 +S 2*^ Z ^C 22 2 ) ) 

T33=2^REAI/C 1 1*G 14C 21 ^2 ) 

T44 = 2-':-RE AL( Cl* ( C 12 *H 1 4C 22 *H 2 ) ) 

T55= -2 *RE A1 ( S 1*C 1 1 *G 1 +S 2*C 21 *G2 ) 

T66 =-2*REAL(C I >< ( S1*C 12 *H 1 +S 2*C 22 2 )) 

CAICULATE STRESSES IN ORTHOTROIIC SHEET 

32 CONTINUE 
S7=XF*DX' T11*D(I)+S7 
S8=XF-< DX>' T22*D ( I+MX)+S8 
S3=XF*DX*T5;-'D(I)+S3 
S4=XF* DX • T44*D ( I+NZ ) + S4 
S5=XF*Da> T55*D( I )+S5 
S6=XF* D::^< T66 *D ( I +I'K ) +S6 

7 OONTIimE 
RETURN 
HID 
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COraSX FUNCTION DBZ(Z,ZO) 

CO MMO N/TO :'/E3 , T 2 , V 1 , SMX , SM Y , G , CO NS , Q , A 1 , SCO K 

THIS FUNCTION IS THE DERIVATIVE OF B(Z,ZO) V/ITH RESISCT 
TO Z 

com EX ZB.CXSR 
CGMELEX Z,ZO ,B1 ,B2,B5,Z0B 
CCNIDEX XIZO ,XI ZOB , XI Z 
ZB=CCN JG(Z) 

ZC3=C0NJG(Z0) 

A=A1 

XIZ=CX3R(Z,A) 

XI ZO=CXSR(ZO,A) 

XI ZOB=CZSR(ZOB,A) 

B1 = (-4'^Q^Z^ ZOB /(Z>^*2-ZOB*^2)**2-4^Z^ (3*ZO*^2*ZOBfZ^*2 
1*ZOB-4>^ZO'^ 3 

2 ) /( ( ZO-Z ) ( ZO+ Z)>^^ 3 ) ) /2 

B2=-((Q^XIZOB* (A-'X2>^ZOB^^'2-2 >^Z'' ^4+A^ '^2*Z*-'^2)/(ZOB>^^4-2 
1*Z^^^2* ZOB 

2*^( 2+ Z - ■>' 4 ) -XI ZO’' ( A"' ^ 2*Z 0^ 2- 2 *Z 4+A * 2* Z >• * 2 ) / ( Z 0 ** 4-2 
1Z-"^2*ZC”^2 

2+Z-'-^’4))/(Z' ^2-A* -2)-+ ((Z0^*2+Z^ Z0-2*k* '2)/(Z0-Z) ^ ZO 
1->'>'2-Z^20 

5-2^A->> ^'2)/(Z0+Z 3+Z^^ (( A-s<-» 2- Z^ Z0)/(Z0-Z)*^2- (Z^ Z0+A^^2 

1)/(Z0+Z)^^2 

4)/(Z"*2-A’ *2)) > (Z0-Z03)/(2*XIZ0))/XIZ 
DBZ=B1+B2 
RETURN 
EICD 

SUBROUTINE RESIGM(Z,S11 ,S22 ,312 ,KTYPE ,STRESSI'i,STRSSSC ) 

COMUTE 3TIGSJES IN ADHEREIIDS DUE TO RH-.OTE STRESSES 

CCM:0 N/TO 1 /E 3 , T 2 , V 1 , SNX , Sl-i Y , G , CO ns , Q , a 1 , scoi: 

DII'Li.'SIOr STRESSi:( 3) ,STRESSC(3) 

CC INI EX Z , GX 3, '' , xN I , or : EGA 3 , DI HI , XX , ZB , SI 
IF(IE’YIE.EQ.1)1 ,2 
1 CONTI :ruE 
ZB=CCNJG(Z) 

:HI =3 T?E 3 -I' ( 2 ) " ( Z/CX SR( Z ,A 1 ) ) /2 . 23 " ( S T.GS 3I.( 2 ) -ST RES 

isr.(D) 

OKEGA. 3=S TRE 3 SI ( 2 ) ' ( ZB / CX SR( ZB,Al))/2.+ .23*(S TEES SN ( 2 ) - 

1STro3S3I ( 1 ) ) 

D1HI=-A1^ ''2/(CXSR(Z,Al)^ ( Z^ -'2-A1>^-2)) 
DlHI=D.HI*3TRESS3(2)/2 . 

SI =a NJG (aCEGAB ) -?II I-Z* Drill 
xSHI=2-REAL(iHl) 

XK=XI'HI+ZB!''DniI+SI 

S1 1 =REAI, ( :C HI- ( 2B*Dl KI+ SI)) 

S22=RE;iJ, ( YJ ) 
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S12=AIMG(Xi:) 

GC’ TO ^ 

2 311=3TRESSGM) 
o22=GT}SSSG(2) 

312=STHES3C!(3) 

3 co:iTii:uE 
RSTURII 
ENG 

COIirLEX FUJfCTION G(Z,v/) 

COKFLEX Z,ZB,V.',V/B 
ZB=00KJG(Z) 

>ffi=00NJG(W) 

G= C LC G ( ',;-Z ) - CL 0G( Z+\fB ) -CL 0 G ( Z+ W) + CLO G ( V/B-Z ) 

RSTUIGI 

CCKREX FUIJCTIOi: H(Z,’.'0 
CCF.TEX Z,ZB,V;,VC ,CI 

ci=aiiLx(o.o,i .0) 

ZB=G0IIJG( Z) 

Vffl=OGIIJG(V/) 

K= CL0G( Z ) + CLO G ( Z+ Vffi ) -CL 0 G ( Z+ VO -CL 0 G ( WB- Z ) 

H=CI>-H 

RETURN 

EI'ID 

SUBRDU TE; E PLA STI C ( SYI ELD , GAD2 , S ,DR ) 

PERFORi: I!:CRL'2SNTA.L PLASTIC ANALYSIS 

CGFM 0 iI/RCT/Z 9 ( 10000) ,D( 100 ) ,NA1F 
COXKGN/ALHES /TA D , GA D 

dileI'jsigi: f( ioo) ,g( ioo),m( 100 ) ,kgp( 100) ,iiyieie( 100) 

I'Z=0 


CGI'ZUTE YIELD STRESS FOR EACH POINT 
CHOOSE CRITICAL ELEi:a'T 


N=2*rAJiF 
NS 0=1'?':: 

SUK=0 

DC 55 1=1, NAIF 
F( I )=0.0 
F( I+NAIF)=0.0 
55 N0x(l)=0 

PCI =99 5999 999 9 9 59. 

12 FCC=9S 999 S 989. 

K}:=:z+i 
DC 11=1, UAL? 

IF(NDi (I).SQ .1)GC TG 1 

ADD LOGIC TO SKIl ALL YIELDED EIEKEn’S 
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A=D(I )^*2+D( I+1'JALP)**2 
B=2*(F(I)^D(I) + P( I+I.IAL?)*D(I+NAI,F)) 

C=? ( I) ''^*2+F ( I+mLF ) X * 2 - SYIilLD^’^2/ 5 . 

TS S T = ( -Bf SQ RT( 2-4 *A* C ) ) / ( 2 . *A ) 

IF (TE3T.lt .FCC)2, 3 

2 K=I 
FCC=TE2T 

3 CONTIIIUE 
1 CONTINUE 

IF(EK.B3 .1 )FCI=FCC 
LL=i:-1 

3UM1=SUM 

3un=sui:+Fcc 

IF(SUI-:.GT.S)5,6 

6 IF(FCC.B2 .999999999. )G0 TO 13 

lT0f'(K) = 1 
NYIELD(IGC)=L 

SAVE STRESS IN EAffl ELEMENT ATYIELDniG 

DC 8 1=1, N 
F(I)=PCC*D(I)+F(I) 

G(I)=FCC-D(I) 

8 CONTIIIUE 

MODIFY EQUATION SET ICR YIELDING OF CRITICAL ELSI-IEilT 


I1=(K-1 )^i:+K 
I4=K+NAI.F 
I2=(I4-1 )*N+I4 

Z9 (1 1 )= Z9(1 1 )+ T AD»^ ( 1 . /GAD2- 1 . /GAD ) 
Z9(I2)=Z9(I2)+TAB^( 1. /GAD2-1./GAD) 

USING AN ITERATIVE METHOD UL'DATE ELEMENT STRESSES 

CAIL GAU5IS(DR) 

IF(I:K.JQ .HATF)GO TO 18 
GO TO 12 

5 IF(n..3Q .1 )15 , 1^ 

13 HUNT 13 
3UM=0. 

13 FC?J-iAT(/>^ THE SOLUTION IS COI'TLETELY ELASTIC^/) 
GO TO 16 

14 CONTINUE 
KK=K-1 

SUM=SUM1 
18 CONTINUE 

tRINT 31,KIC,aJM 

31 FORMAT( /no* ELEI'CHNTS HAVE YIELDED AT ^E10.3/) 

16 CONTINUE 
1001 F0RMAT(8E10.3) 
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DO 26 1=1, mF 
D(I)=(S-SUM)-D(I)+P(I) 

D( I +NALF) = ( S-SUM )*D( I+NALP ) +F( I+NALF ) 

26 CONTINUE 
IF(KK.GT .1 )25,29 

25 HUNT 27 

27 I0RI:AT(/* YISTiDSD EIEI'.ENTS */) 

HUNT 28, (NYIEiD(l),I=1,KK) 

28 FC'RMAT( 12I10 ) 

29 OONTIIIUE 
IRINT 30, PCI 

50 PCRHAT(^ THE YIELD MACROSOOPIC STRESS IS *E11.3/) 
RETUI^J 
END 

SUBRDUTKE GAUSS (DD) 

GAUSS SEIDEL METUOD FOR SOLVING HICREMUITAL 
n-ASTTC SCLUTiai 

DIMENSION DD( 100 ) , ASAVS( 100 ) 

COKI*DN/ADHES/TAD ,GAD 
C0MM01I/R0T/Z9( 10000) ,D( 100 ) ,NALF 
EPS=.005 
N=2*NAIF 
NSQ=mN 
ITMAX=20 
DO 33 1=1 ,N 
K=(I-1 )*N+I 
ASTAR=Z9(K) 

ASAVE( I )=ASTAR 
DO 3 J=1,N 
II=(J-1 )^N+I 
Z9(II )=Z9(II )/ASTAR 
CON TUTU E 

DD(I)=DD(I)/ASTAR 
33 CONTINUE 

DO 9 ITER=1 ,ITI-LAX 
KFLAG= 1 
DO 7 1=1, N 
XSTAR=D(I ) 

D(I)=DD(I) 

DO 3 J=1 ,N 
II=(J-1 >»N+I 
IF (I .35 . J) GO TO 5 
D(I)=D( I)-Z9(II)^D(J) 

5 CONTINUE 

IF( ABS(XSTAR-D(I ) ).LE.EPS)GO TO 7 
KPLAG=0 
7 CONTINUE 

IF(KFLAG.NE.1 )G0 TO 9 
GO TO 1 
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9 CONTINUE 
PRINT 204 
1 GONTH'IUS 
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RECOIBTRUGT Z9 AND DD MATRIX K)R USE ON PUOUTRE 
ELEMHITS THAT HELD 

DO 350 1=1, N 
K=(I-1 + l 

ASIAR=ASAVE( I ) 

DO 50 J=1,N 
n=(J-1 >11+1 
Z9(II )=Z9(II )^ASTAR 
50 CONTINUE 

DD(I)=DD(I)*AS0AR 
330 CONTINUE 

204 FORIAT(/* CAUTION THE GAUSS SEIDEL DID NOT CONVERGE^--/) 
RETURN 
HID 

FUNCTION YYX(X,B,A) 

COMMON /DOND/F , n , P2 , XKF , XKUN 3 , XKS T IFF 
IF(X.B3.0)G0 TO 4 
IF(X/A.(2:.1 )1,2 

2 YYX=I^( (1.-(X/A)^ *^Pl)**(l./r2)) 

GO TO 5 

1 YYX=0. 

GO TO 5 

4 YYX=B 

3 CONTINUE 
RETURN 

FUIICT ION YYD ( X , B , A , P) 

IF(X/A.GE.1 )1,2 

2 YYD = I^((1.-(X/A)’ )-*‘-(l./P )) 

GO TO 5 

1 YYD=0. 

3 CONTINUE 
RETURN 
END 
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